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1£. ABSTRACT 

A computer mettiod is presented, using the "inside-out" approach, for predicting 
the structure of the disturbed zone near a moving bo<fy in space. The approadi uses 
fewer simplifying assumptions than otiier available methods, and is applicable to 
large ranges of the values of bo<fy and plasma parameters. Two major advances 
concerning 3-dimensional bodies are that (a) tiiermal mc^ions of ions as well as of electrons 
are treated realistically by foUowii^ their trajectories in ttie electric field, and (b) the 
technique for achieving self-consistency is promising for very large bodies. 

Three sample solutions are obtained for a disk-shaped bo<fy, charged negatively to a 
potential 4kT/e. With ion Mach number 4, and equal ion and electron temperatures, 
the wakes of a relatively smaU bocfy ( radius 5 Debye lengths) and a relatively large 
bo<fy (radius 100 Debye lengths) both begin to fill up between 2 and 3 body radii downstream. 
For the large body there is in addition a potential well (about 6kT/c deep) behind the body, 
hicreasii^ the ion Mach number to 8 for the large boch'^ causes the potential well to become 
wider and longer but not deeper. For the large body, the quasineutrality assumption is 
validated outside of a cone-shaped region in the very near w'ake. For the large as well as 
the small body, the disturbed zone behind the body extends transversely no more than 
2 or 3 bodj' radii, a result of significance for the design of spacecraft boom instrumentation. 
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FOREWORD 


This report presents the results of a study to access the plasma 
disturbance created by a large body of Space SSiuttle dimensions 
traversing the ionosphere. This study was conunissioned by the 
Plasma Flow and Interaction Section of the Atmospheric Magnetospheric 
and Plasmas in Space (AMI^) science deHnition working group. The 
study was performed by Lee Parker, Inc. of Concord, Mass, for the 
NASA, Marshall Space Flight Center under the direction of W. R. Roberts 
of the AMPS Task Team. 
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CHAPTER 1 


IHTRODUCTION 

The problem of theoretically calculating the structure of the dis- 
turbed plasma {Frequently referring to the wake and/or sheath) around e 
moving body in space is equivalent to that of solving a complicated sys- 
tem of coupled nonlinear partial differential/integral equations. The 
equations consist of the Vlasov (collisionless Boltzmann) equations for 
the ions and electrons, and the Poisson equation relating the electric 
field to the distributions of ions and electrons. The difficulty is 
essentially a numerical one because analytic solutions are not possible 
(for cases of interest), and there is no unique approach. In cases of 
stationary bodies (Parker, 1973 and 1975), as well as moving bodies 
(other theoretical references of this report), combinations of numerical 
techniques (finite differences, iteration, quadratures, etc.) are required 
for treating various parts of the problem. For either stationary or mov- 
ing bodies, the choices of techniques and their use to achieve consistent 
solutions for any given set of physical parameters (defining body and 
plasma) have never been obvious. Innovations are frequently required. 

The purpose of this report is to review some of the available techniques 
for a moving body (with emphasis on the wake), to describe in detail a 
new combination of techniques which appear to be reasonably successful 
over a large range of the physical parameters, and to present sample 
solutions as well as the implementing computer program. 

Various approaches which have been used for this type of problem are 
summarized in Chapter 2. In all such calculations simplifying assumptions 
are made. The customary ones are: 

Collisions negligible. 

Geomagnetic field negligible. 

Simple geometry (sphere, disk, cylinder, etc.) 

Simple surface reactions (usually, charged particles are 

neutralized). 

Prescribed surface emission (usually none, but simplified 



photoelectron and secondary-electron emission are 
includable). 

Conducting body (usually perfectly conducting, but 
finite conductivities are includable). 

Steady state. 

These assumptions may be questioned (for example the neglect of time- 
dependent phencmiena), but they may be at least partially relaxed by employ- 
ing known techniques to generalize the calculations. In the interest of 
achieving reasonably economical calculations witnin the limits of available 
computers, the above assumptions in their usual form are adopted in the pre- 
sent work. 

The techniques and computer program described in Chapters 2 and 3 and 
in the Appendices have been developed to solve the coupled Poisson-Vlasov 
system of equations to obtain distributions of ion and electron density, 
and potential, about 3-dimensional bodies (with axial sytrenetry about the 
direction of plasma flow). The program uses the "inside-out" method devel- 
oped by the author in 1%4, which follows ion and electron trajectories 
backward in time from the point in space at which it is desired to know the 
velocity distribution, to their origin in the undisturbed plasma where the 
distribution is known. 

Briefly, the present approach (see Chapter 2) differs from that of 
Call (1969) and Martin (1974) by including both the ion and the electron 
thermal motions, whereas Call and Martin represent the distribution of ions 
by a cold beam, and of electrons by the Boltzmann factor. The approach 
differs from that of Taylor (1967) in that (a) it is applied to 3-dimensional 
bodies whereas Taylor treats an infinitely-long cylinder, and (b) the Poisson 
and Vlasov calculations are cycled until self-consistency is achieved, where- 
as Taylor's calculation is terminated after the first cycle. The approach 
differs from that of Grabowski and Fischer (1975) because they (a) assume 
that quasineutrality holds throughout space (which is invalid in the very 
near wake), and (b) apply their method to an infinitely-long cylinder. Dif- 
ferences with other methods are outlined in Chapter 2. The most similar 
calculation previously done was for an infinitely-long cylinder by Fournier 
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(1971), using the inside-out niethod. In general, the present approach uses 
fewer simplifying assumptions and is thus applicable to a larger range of 
parameters than other available methods. 

Two major advances are represented by the present program, as opposed 
to previous approaches, particularly with regard to wakes of 3-diraensional 
bodies: 

(1) Thermal motions of ions as well as of electrons are treated 
realistically by following their trajectories in the electric field. 

(2) The technique for achieving self-consistency is promising for large 
bodies many orders-of-magni tude larger than the Debye length (the Shuttle- 
Orbiter or the moon, for example). 

Solutions may be obtained with reasonable amounts of computer time by 
judicious choices of grid points and other numerical parameters. The 
method can be extended to include an arbitrarily-shaped body (presently a 
body of revolution). 

The structure of this report is as follows. 

Chapter 2 comprises a review and summary of previous approaches, clas- 
sified on the basis of how they treat the Vlasov problem (calculation of ion 
and electron densities and currents). In particular, the inside-out method 
is treated in detail. The computational method for number and current den- 
sity quadratures is given in Appendix A. (Throughout the report the words 
"orbit" and "trajectory" are used interchangeably.) 

In Chapter 3, the method of self-consistent solution by Poi sson-Vlasov 
iteration is treated. The method of solution of the Poisson problem by fin- 
ite differences is described in Appendix B. The "ion-density option" dis- 
cussed in Chapter 3 is appropriate for the large-body problem (see also 
Appendix B and Chapter 4). 

The FORTRAN listing for the computer program, and the descriptions of 
input and output data, are given in Appendix C. 

Chapter 4 presents numerical solutions for three sample disk problems, 
showing the effects of changes in body size and in ion Mach number. The 
results are presented in the form of transverse profiles of ion c'^nsity, 
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electron density, and potential, in the wake. (Data on the sheath in 
front of the disk are available but are not given in this report.) The 
key results are the following. 

For a relatively small body (radius = 5 Debye lengths), at potential = 

-4 kT/e and ion Mach number = 4, there is no prominent wake structure such as 
large-amplitude "bumps" (enhancements/depletions) in the ion or electron density. 
The wake becomes filled in between 2 and 3 body radii downstream, and there 
is no potential well. The wake disturbance is essentially confined to a 
region of length in the axial dimension about 4 radii, and transverse 
radial dimension about 1.5 radii. 

For a large body (radius = 100 Debye lengths, i.e., larger than has 
been previously treated realistically), and for the same body potential and 
ion Mach number as above, the wake begins to fill up again between 2 and 3 
body radii downstream. The wake disturbance extends more than 6 body radii 
downstream, but transversely only between 2 and 3 radii. There is a poten- 
tial well near the wake surface of the disk, and quasi neutrality is valid 
outside of a cone-shaped region near the wake surface (with the disk forming 
the ba»e of the cone). 

For the same large body, but with ion Mach number 8 instead of 4, the 
dimensions of the wake-disturbance region are not significantly changed, 
but the filling-up occurs further downstream. The potential well becomes 
wider and longer, although the depth is similar. In addition, there seems 
to be a central core of essentially ambient density along the axis, for 
both ions and electrons. 

To more comprehensively establish the practical applicability of the 
present computer method to future AHPS/Spacelab missions, it would be of 
interest to compare theory and experiment for cases where in- si tu and labora- 
tory simulation data are available. At present there are more laboratory 
results (Oran et^., 1975; Fournier and Pigache, 1975), than in- «;itu results 
(Henderson and Samir, 1967; Sa~.ir et ^. , 1973). However, it is presently 
still difficult to simulate ion transverse velocity distributions in the 
laboratory, and the effective ion temperature is frequently low Since vari- 
ous ratios of ion temperature to electron temperature may be treated by the 
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program with relatively minor modifications, computations should be made 
with "cool" ions to facilitate lab- theory comparisons. The present theo- 
retical nwdel should also be compared, using selected ionosphere-magnet- 
osphere example problems, with other available theoretical models. The 
other models may be less realistic but they may have advantages of rela- 
tive simplicity and economy; they may also be "calibrated" through such 
comparisons. 

The present computer method gives information regarding the dimensions 
of the disturbed zone about a body. Information of this kind should be 
useful for estimating the lengths of booms to be deployed, for example, 
on the Spacelab to keep outboard instrumentation outside the disturbance 
created by various structures on the spacecraft. In this sense, the com- 
putations may be regarded as a phase of a feasibility study. 
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CHAPTER 2 


APPROACHES: REVIEW AND SUMMARY 

All of the approaches to the body-in-a-plasma problem have in common 
the following elements. The quantities to be computed include (a) the 
potential distribution and (b) the ion and electron density distributions. 
One may also include the associated surface current densities. The equa- 
tions to be solved simultaneously are (a) the Vlasov equation for ions, 

(b) the Vlasov equation for electrons, and (c) the Poisson equation. The 
solutions of the Vlasov equations (velocity-distribution functions) are 
used to compute number densities (and surface current densities). The 
number-density distributions become input to the (right-hand side of the) 
Poisson equation v^hich yields the potential distribution. Finally, an 
iterative procedure is used for self-consistency, v/herein the density and 
potential distributions are successively cycled until satisfactory conver- 
gence has been achieved. 

The steady-state Vlasov equations for ions and electrons state that 
the velocity-distribution functions remain constant along particle tra- 
jectories. With the electric field assumed given (numerically in terms 
of a spatial grid about the body), solving the Vlasov equations means 
formally that one determines, from the shapes of the trajectories, the 
ion and electron velocity distributions at the grid points. The trajec- 
tories relate local velocities at a given grid point to those at infinity. 
Through these relationships, the ion or electron number density at the 
point may be evaluated by a velocity- integral over the local velocity 
distribution. Similarly, the current density may be evaluated at desired 
locations (usually the body surface). 

It is convenient to classify the various theoretical approaches on 
the basis of hew they treat the trajectory part of the Vlasov problem. 
"Inside-out" methods follow the trajectories backward in time into the 
undisturbed plasma, while "outside-in" methods follow the trajectories 
forward, in the direction of physical motion of the particles. (In an 
outsidc-in method, the velocity-distribution function is not calculated; 
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rather, the density is evaluated directly.) "Other" methods denote approxi- 
mations where trajectories are not followed at all. The three approaches 
are discussed below. 

There exists as yet no systematic comparison of the results of the 
various approaches with one another. 

Before discussing the various approaches, we may define here the 
parameters of interest: 

Plasma Parameters 

n^ = unperturbed number density at infinity 

T^,Tg = ion, electron temperatures 

m- = ion mass (electron mass not required) 

Aq = electron Debye length 

Body Parameters 

r„ = characteristic dimension 
0 

Vq = velocity 
= body potential 

4>o = " dimensionless body potential 

M = V /m./2kT. = ion Mach number (electron Mach number assumed 

negligible) 

Xq = = Debye number 

Hencefortii aVi lengths are to be considered normalized by r^. Thus, Xq 
will denote the dimensionless Debye number. Potentials are normalized by 
kT /e, so that <})(r) denotes the dimensionless potential at the spatial 
point r. Number densities are normalized by n^, so that n(r) denotes the 
dimensionless density at r. In the calculations involving integrations 
over velocities, v will denote a velocity normalized by the value of 
/ 2kT/m associated v/ith the particles of interest. Similarly, E will 
denote total energy normalized by kT, Velocity-distribution functions 
(denoted by f) will be normalized by n^. For a given body geometry, there 
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are four dimensionless physical parameters of interest, namely, x^, 
and M, along with the temperature ratio T./T^. Table 2-1 shows a sampling 
of the parameters used in various previous calculations. 


2A. IHSIDE-OUT METHOD 

Consider a single species of (charged) particle, i.e., ions or elec- 
trons, The electric field is assumed to be known. In order to compute 
the number density n(r) at the point r, one must evaluate the triple 
integral over velocity space: 


n{?) 



) dv^dVydv 


z 


(2-1) 


where f(r, v) is the distribution function which satisfies the Boltzmann 

equation for the given species of particle, f is the radius vector of the 

space point of interest, and v is the local velocity of a particle at r. 

The velocity-volume element is written as if cartesian coordinates were 

being used, but the product dv dv dv is intended to symbolize an arbitrary 

X y z 

coordinate system. Similarly, in order to compute the collected current 
density at points on the surface of a body, one must evaluate at each point 
a triple integral over velocity spa'.e of the form 


j 




dv dv dv 
n X y 


( 2 - 2 ) 


where v^ is 
n 

at the point 


the component of the particle velocity normal to the surface 


r. 


The problem is thus to evaluate f. Since the problems of interest 
are assumed to be collisionless and constant in time, the distribution 
function f satisfies the steady-state Vlasov (or collisionless Boltzmann) 
equation, namely. 


V • vf + a * Vyf = 0 


(2-3) 
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TABLE 2-1 


PARAMETERS ADOPTED IN PREVIOUS WAKE CALCULATIONS 



Mach Number 

Debye Number 

Ion-Attractive 
Body Potential 

Fournier (1971) 

1. 6, 10 

1 2 
10 * 3 

-3, 0. 1. 2.75, 
6, 40 

Davis and 
Harris (1961) 

6 -»■ 7 

1 1 
25 * 10 

20 ^ 1000 

Call (1969) 

1 8 


0 ^ 40 

McDonald and 
Smetana (1969) 

0 ->■ 3 

1, si . 10 

10, 25 

Maslennikov and 
Si gov (1964) 

2. 7 

1. 12 

0 -»■ 40 

Liu and Jew 
(1968, 1969) 

4, 8 


1. 5 

Kiel et al. 
(1968) 

5, 8 

1 1 1 
TO’ T0O’ 10 

3 

Martin (1974) 

4 ^ 10 

1 

5, 9 

Grabowski and 
Fischer (1975) 

0 1.4 

0 (quasineutral) 

irrelevant 

Taylor (1967) 
(first order only) 

6 

2/3 

0.25 -V 20 
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where a is the vector acceleration of a particle passing with velocity v 
through the point r. The gradient operators 7 and 7 operate on the com- 
ponents of r and of v, respectively. Equation (2-3) states that f is con- 
stant along a particle orbit, which is characterized by the constants of 
the motion. In a general electrostatic field (here assumed given) whose 
sources are volume and surface charges, the total energy E is conserved, 
where the dimensionless E is defined by 

E - + Hr) (2-4) 

and <^(r) is the dimensionless potential energy of the particle at r. 

With <J>(r) a known function of f, one may evaluate the integrals in 
Eqs. (2-1) and (2-2) by following orbit5. backward in time with trajectory 
calculations to a point where f is known. For example, in the case of a 
body immersed in a plasma, f is assumed to be known at infinity (where (». 
vanishes), and is assumed to have at infinity a prescribed energy distri- 
bution, such as a Maxwellian with drift, or a more general distribution. 
Also, f is assumed to be known on the surfaces of electrodes. If a sur- 
face emits particles, its distribution function must be prescribed. If 
the surface absorbs without re-emitting charged particles, the distribu- 
tion function (of emitted particles) is prescribed to be zero. Thus, f 
is discontinuous in velocity space: That is, the physically-possible 
velocity space (at the point r) is divided into two domains, namely, the 
domain of orbits which have come to r from infinity, and the domain of 
orbits which have come to r from electrode surfaces. In the latter 
domain, f vanishes if there is no emission. There rore, f is discontinu- 
ous on the boundary between the two domains in velocity space. The shape 
of the boundary between the two domains depends, of course, on the geom- 
etry and the potential function in and it is the heart of the problem 
(a) to determine the boundary of the domain of orbits coming From infinity, 
and (b) to evaluate the integrals Eqs. (2-1) and (2-2) over that domain of 
velocity space. 

In practice, one need not in general determine explicitly the boun- 
dary of the domain in velocity space of orbits coming from infinity. 
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Rather, one may follow a large number of orbits backward in tin« (computa- 
tionally), and evaluate the moment integrals, Eqs. (2-1) and (2-2), auto- 
matically from the results of the orbit-following, it may, hoviever, under 
scMne ci rcimistances be more accurate and efficient to determine this boun- 
dary. To do so would complicate the canputer prograRining. 

For a Haxwell’an distribution driftirj with Mach number H, the dimen- 
sionless velocity-distribution function at infinity may be written: 


f 


, -(v ^ - 2Mv ) 

I • «o 

-^372 " 


O 0 

^ -(^+v +M 

^372 ® 




(2-S) 


(velocities in units of »'?Fr/m , 

= axial component of velocity) 


2 2 

v^here v = v + *» may be identified with the total energy E, and v, 

00 

with times the cosine of the angle between v and the axis. The 

<o 

moment integral (2-1) for number density may be approximated by a quad- 
rature sum as follows: 


n 



d v 


■I I I *Uk «ijk 

i j K 


(2-6) 


where d v is a short-hand notation for the element dv dv.dv,, and 6 is a 

X jf z 

cut-off (or step) function, equal to unity or zero according as the tra- 
jectory is found to come from infinity or the body surface, respectively. 

In the sum, the three indices refer to discrete values of three components 
of velocity, where the values are chosen in accord wich a quadrature scheme 
(Gaussian), and the coefficients A... are proportional to the associated 
weights. Each term in the sum represents an individual trajectory. A sim- 
ilar sum is obtained for the current density (see Appendix A). 

Figure 2-1 indicates schematically how one of the trajectories (with 
indices i, j, k) from the sum in Eg. (2-6) is traced backward from the 
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To grid boundary (infinity) 65^,^ = I , and 

use computed v 









— local velocity (v 


1 

1 

p 







'ijk 


n is evaluated 


To body surface = 0 



Evaluation of 


for (i,j,k)-th trajectory by following (reversible) 


trajectories backward in time. 


FIG. 2-1. INSIDE-OUT METHOD 
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point P (usually, a grid point), and found either to reach the body sur- 
face, or to reach “infinity" at the boundary of the grid. 

This constitutes the inside-out im:thod of solution of the Vlasov 
problem. Further details including the discrete velocities and coeffi- 
cients of the sum are given in Appendix A. 

The advantages and disadvantages of the n:ethod are: 

Advantages of Inside-Out Method 

1. Density points can be chosen individually and at random. 

Hence the nethod is flexible. 

2. Suitable for electrons as v#ell as ions. 

Disadvantage (relatively minor) 

Information carried by trajectories is lost upon moving to 

another density point. Hence the calculation tends to be 

time-consuming. 

The inside-out n«thod was developed by Parker (1964), and has subse- 
quently been used by Fournier (1971) and by Grabowski and Fischer (1975) 
to calculate the v/ake of an infinitely-long moving circular cylinder. 
Grabowski and Fischer also assumed complete quasineutrality, thus restrict- 
ing the generality of their method. It was also used by Taylor (1967) for 
the wake of an infinitely long cylinder of rectangular cross-section (a 
"thick strip"), but the calculation was not carried beyond the first itera- 
tion, and is therefore not self-consistent. Parker and Whipple (1967, 1970) 
have used the method for two-electrode probes on a satellite, and Parker 
(1970, 1973) has used the method for two-electrode rocket-borne and labora- 
tory probe systems, and for the problem of a small probe in the sheath of 
a large electrode. 


2B. OUTSlDE-IN METHODS 

Outside-in methods may be divided into two types, the method of flux- 
tubes and the method of superparticles or weighted deposition. The two 
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types of outside-in methods are illustrated in Fig. 2-2. The trajectories 
are injected from the outer boundary of a grid. 

In the method of flux tubes, the flux of particles in a tube is con- 
stant . The tube is defined by two neighboring trajectories. Since the 
cross-sectionai area of the tube is known from the trajectory calculation, 
and the particle velocity is also known, the particle density may be deter- 
mined at any point in the tube (as indicated in the figure). The density 
is usually assigned to the nearest grid point along the path of the tube. 

The advantages and disadvantages of the flux- tube outside- in method 

are: 

Advantages of Method 

A relatively fast calculation. 

Disadvantages 

1. Invalid if trajectories cross or reverse direction. 

Hence the region near the body's wake surface cannot 
be treated. 

2. Suitable only for an axisyNinetric body with cold ions 
in a beam. 

The flux-tube technique has been used by Davis and Harris (1961) for 
the cold-ion wake of a sphere, by Call (1969) for the cold-ion wakes of a 
strip, disk, infinitely-long cylinder, and sphere, by Martin (1974) for 
the cold-ion wakes of a strip and a disk, and by McDonald and Smetana (1969) 
for the wake of an infinitely-long cylinder in a monoenergetic-ion plasma 
with drift. 

In the method of weighted deposition, the space is divided into cells, 
with each cell associated with one of the grid points. The contribution of 
a trajectory to the density in the cell is proportional to the time spent 
in passing through the cell. 

The advantages and disadvantages of the method of weighted deposition 

are: 
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METHOD OF FLUX TUBES 
Tube defined by two neighboring trajectories 



N = number entering 
tube per second 
at outer boundary 
of grid. 

= nAv = constant 


A, A' , A" = cross-sectional area at points P, P*, P" 

V, v', v" = local speed at P, P' , P" 

n n r, - N 

"p ■ Av * "p* “ A'v' ’ V " 

(Method assunres n = constant on cross-section of tube) 


METHOD OF SUPERPARTICLES 
OR WEIGHTED DEPOSITION 



Contribute to 
density in cell 


Contribution to density in cell proportional to time (At = As/v) spent in cell. 

*n = N At ^ N AS 

~ (volume of cellj v x (volunx; of cell) 

(method assumes n = constant within cell) 

FIG. 2-2. OUTSIDE-IN METHODS 
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Mvantaqes of >tethod 

1. N*' difficulty with trajectory crossings. 

2. Related to and adaptable to t1me>dependent computer 
simulation. 

Disadvantage 

Many trajectories needed tor good statistics within cells. 

This method was studied by Parker (1964) for a monoencrgetic-ion dis- 
tribution with drift* and was used by Maslennikov and Sigov (1965) for the 
cold-ion wake of a sphere. 

2C. OTHER METHODS 

Other methods include approximate treatments which avoid trajectory 
calculations. Liu and Jew (1%8* 1969) assumed that the ion axial compo- 
nent of velocity is constant. They then determined limiting trajectories 
for the density integral by further approximations, namely, an additional 
assumed approximate constant of the motion, evaluated using the local field 
in the vicinity of the point in question. They applied their method to the 
wakes of a sphere and a cylinder. Kiel, Gey, and Gustafson (1968) treated 
the wake of a sphere, assuming neutral ion trajectories (straight-line paths, 
neglecting the electric field). They also assumed that the electron densi- 
ties were given by approximate formulas designed to include the effect of 
the potential barrier in the wake. For the wake of a strip, a disk, and an 
infinitely-long cylinder, Gurevich ^ aT. (1969) assumed quasineutrality, 
with ion and electron densities both equated to the Boltzmann factor. In 
addition, they assumed that the ion axial component of velocity is constant 
and that the ion Mach number is large. 
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CHAPTER 3 


SELF-CONSISTENT SOLUTIONS BY POISSON -VLASOV ITERATION 

The "inside-out" method for obtaining ion and electron densities* in 
a given electric field defined by the values of the electric potential at 
a chosen set of grid points, has been described in Chapter 2, with compu- 
tational details given in Appendix A. This constitutes the "Vlasov prob- 
lem." The Vlasov problem must be solved separately for the electrons and 
each species of ions (when there is more than one). In going from one 
species to another, or to electrons, the potentials are multiplied by the 
appropriate scale factor. 

How the electric field is obtained when the ion and electron densi- 
ties are given is discussed in detail in Appendix B. Here, the Poisson 
equation is replaced by a set of difference equations based on the chosen 
set of grid points, with one equation for each unknown potential at a grid 
point. The derivation of the coefficients of the unknown potentials in 
the difference equations, and the method of solution, are given in Appendix 
B. The system of simultaneous equations for the unknov/n potentials is 
solved by a relaxation procedure. This constitutes the "Poisson probl«n." 

The boundary conditions for the potentials in the Poisson problem 
are as follows. At points representing the body surface, the normalized 
potential is fixed at the chosen value 41^. At the external (boundary) 
points of the grid, where "infinity" is represented on the computer, a 
"floating" condition is optionally used, namely, a linear relation between 
^ and 34i/3n, the normal component of V(|>. The exact relation of 4» to 3(ji/3n 
is not important when the external boundary of the grid is sufficiently far 
away. (For the calculations to be repoi^ted, the assumed relation was the 
same as for a Coulomb potential.) In any case, either the fixed condition 
^ = 0 or the floating condition will give the same results, provided the 
grid boundary is moved sufficiently far out. The effects of various types 
of boundary conditions representing "infinity" have been studied by Taylor 
(1967), and by Parker and Sullivan (1974). In general, the floating con- 
dition appears to be computationally more efficient than the fixed one. Of 
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course, the floating condition becomes ideal when the true relation 
between 4 and a^/an is used, but this requires that the asymptotic form 
of the solution be known in advance. (See, for example, Parker and Whip- 
ple (1970).) 

An Iteration method is used for computing self-consistent charged- 
particle and potential distributions. This is herein referred to as the 
"Poisson-Vlasov iteration." Two principal options are employed for this 
procedure. In one of the options, the "charge-density" option, the space 
charge is initially and arbitrarily assumed to be zero. For this case, 
one obtains the Laplace (space-charge-less) electric field from the Poisson 
problon. This is the "zero-order" potential distribution, which becomes 
input to the Vlasov problen. The resulting solution of the Vlasov probl^ 
yields the ion and electron densities at the grid points, which ai^ com- 
bined to make "zero-order" charge densities. These become input to the 
next Poisson problem, which then yields the "first-order" potentials, and 
so on. In this procedure one usually "mixes" successive charge-density 
iterates to in^rove stability. Otherwise the process can "blow up." One 
can also mix potential iterates rather than densities if desired. The 
dependence of the stability and convergence of the above procedure on the 
mixing parameter have been studied analytically by Parker (1970) and 
Parker and Sullivan (1974). (No other analysis of this type has been pub- 
lished to the author's knowledge.) This (charge-density) option is most 
effective when the spatial region of interest is not many Debye lengths 
across. The analysis shows that one can (probably always) choose a mixing 
parameter sufficiently small to ensure convergence, but at the expense of 
additional iterations. 

In the other option, the "ion-density option," the ion density distri- 
bution alone is assumed initially. Initial guesses which can be employed 
include (a) zero ion density everywhere, (b) unit ion density (the ambient 
value) everywhere, and (c) the neutral ion density which obtains when there 
are no forces. Whichever choice is made for the initial guess is desig- 
nated the "zero-order" ion density. Now if one assumes the electron den- 
sity to be given by the Boltzmann factor exp(4), the Poisson equation may 
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solved, holding the ion densities fixed, but regarding both the potentials 
and the electron densities at the grid points as unknowns. This is a non- 
linear problem, which is solved by a modification of the relaxation pro- 
cedure used for the "charge-density" option. The new procedure is an 
important advance since the iteration is not as sensitive (tending to blow 
up) to small Debye numbers as in the charge-density option. Thus, very 
large bodies (in multiples of the Debye length) can be treated. This has 
been the method used to obtain the results reported In Chapter 4. Similar 
ideas have been used by Call (1969) and Fournier (1971), but these workers 
have not treated large bodies. 

The assumption that the electron density is given by the Boltzmann 
factor becomes invalid when the body surface potential is near zero, or 
when there is a potential barrier in the wake such chat the wake potentials 
are more negative than the surface potential (causing electrons to be 
attracted to the surface rather than repelled from it). In this case it is 
still possible to use the ion-density option, with its large-body capability, 
provided that, within each cycle, where the ion densities are held fixed, a 
"minor iteration" is carried out such that the electron densities are com- 
puted by trajectory calculations. 
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CHAPTER 4 


SAMPLE RESULTS 

Calculations were matte for three sain)le problems, using the computer 
program listed in Appendix C, based on the theory of this report. The 
results presented here are preliminary in the sense that they are intended 
as an Illustration of the capability of the program, rather than a system- 
atic study. The body is assumed to be a circular disk with its plane nor- 
mal to the flow, and the problem is specified by the dimensionless physical 
paraneters M, and Xq, defined by (Chapter 2): 

♦o = 

M = mvQ^/2kT 

Xq = ^/r^ 

where T is either the ion or electron temperature (assuming equal temper- 
atures), is the disk potential, is the disk velocity, m is the ion 
mass (M is assuned to be negligible for the electrons), r^ is the disk 
radius, and ^ is the dimensional Debye length. 

Numerical parameters for the calculations include 89 grid points, 
distributed mostly in the wake region, and in most cases 512 trajectories 
per grid point (8 values each for the polar and azimuthal angles, and 8 
values for the energy; see Appendix C). 

The potential was set to zero on the downstream computational boun- 
dary, and was allowed to "float" on the upstream and side boundaries. 

(The boundary conditions at the various outer grid surfaces can be either 
fixed or floating.) The downstream boundary for the present calculations 
was set at 6 radii, i.e., beyond the Mach number of radii for the two prob- 
lems with M=4. The electron density n^ was assumed to be given by the 
Boltzmann factor exp(<j»). This is reasonable for itip^-4 on the surface, and 
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leads of course to computer economy (by avoiding trajectory calculations 
for electrons). However, it must be emphasized that this does not repre- 
sent an essential restriction; the program is specifically designed to 
compute rtg, as well as n^. , realistically when necessary, by trajectory 
calculations. Moreover, in cases where a potential v<ell occurs in the 
wake near the surface, the Bo Itzmann-f actor assumption becomes invalid 
and the trajectories must then be computed, at least for points near the 
surface. 

With the Boltzmann-factor assumption, an option is available in the 
Poisson-solution part of the program. This is the "ion-density" option 
(Chapter 3) which includes the Boltzmann factor in the equations for the 
potential distribution; the equations thereby become nonlinear rather 
than linear. This is a simplified case of a possible general technique 
where, during each "major" iteration cycle in which the ion densities 
are held fixed, self-consistent potentials and electron densities are 
simultaneously determined. This technique is as yet in an experimental 
stage, but it seems promising in that it may produce solutions with 
sonable costs for large-body problems; in such problems, the conventional 
Poisson-Vlasov iteration based on the "charge-density" option (Chapter 3) 
becomes expensive (Parker and Sullivan, 1974). A disadvantage of the ion 
density option, hov/ever, is that its convergence properties are not 
understood; therefore its costs are difficult to predict. This is in con 
trast to the case of the charge-density option where an analysis is avail 
able (Parker and Sullivan, 1974). 

The three calculations to be described next were all made with the 
ion-density option. The cases are: 


(a) 

♦o = 

M = 4, 

Xq = 1/5 

(b) 

♦o = 

M = 4, 

Xjj = 1/100 

(c) 

♦o = 

M = 8, 

Xp = 1/100 


Transverse profiles of normalized ion density (n.), electron density 
(n^), and potential (<|>) in the wake region downstream of the disk are 
shown below for the three cases. The profiles are in transverse planes 
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at various distances downstream, and all lengths are normalized by the 
disk radius, (z denotes axial distance downstream, In radii, with z=0 
defined as the plane of the disk, and r denotes transverse, I.e., radial, 
distance from the axis.) The profiles are presented In Figs. 4-1 to 4-3, 
and In Tables 4-1 to 4-9. 
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4.1 Case = -4, M = 4, and Xq = 1/5 

For this case (Fig. 4-1, and Tables 4-1 to 4-3), the parameter values 
illustrate a c"* -s of problems of physical interest, such as a small TAD 
at high altitudes, or a probe mounted on or near a spacecraft. 

Twelve major iterations (Poisson-Vlasov cycles) were computed, in 
which successive ion-density iterates were mixed, with a mixing parameter 
0.5, starting with uniform ambient ion density as an initial guess. In 
the last four iterations, 8192 trajectories were used at points at and 
near the wake axis. This increase of the number of trajectories by a 
factor 16 was made to increase the accuracy in the investigation of pos- 
sible structure (enhancements or depletions in ion or electron density) 
near the axis. 

Figure 4-1 shows three sets of profiles, one for n^, one for n^, 
and one for WHiiin each set, the profiles are arranged vertically in 
order of increasing axial distance z. There are eight values of z, namely, 
z = 0.2, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0 and 6.0. Each profile is constructed 
using nine values of r, namely, r = 0, 0.1, 0.3, 0.6, 0.8, 1.0, 1.2, 1.5, 
and 2.0, with straight-line segments connecting the values of the func- 
tions (n^ , n^, or <|>) computed at these points. (The nine values of r and 
eight values of z are the coordinates of the 72 grid points chosen to 
reprccant the wake region of this problem, for which a total of 89 points 
w. used.) 

The <j)-profi1es (right side of figure) and the n^-profiles (middle 
of figure) are the. 12Lb" order iteration values. These <|)-values are also 
given in Table 4-:j, and the n^-values are given, in parentheses, in Table 
4-1. On the left side of the figure, there are two sets of n, -profiles, 
one labelled "A", and the other unlabelled. The "A”-profi les for n. are 
the llth-order values (Table 4-2), from which the 4>-profiles (and n^- 
profiles) in the figure are derived. The unlabelip-' profiles fo*^ n. are 
the 12th-order n. -values (Table 4-1) resulting from traiectory calculo- 
tions using the i}>-profiles (Table 4-'-}. The juxtaposition of the "A" and 
unlabelled n^-profiles iridU-'. ' • extent to which the Poisson-Vlasov 
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Iteration has converged. At z = 1 and below, the two profiles are suf- 
ficiently close to be considered well converged for the present purposes. 
Tables 4-1 and 4-2 give the numerical values of the 12th and 11th orders 
of n^ (unlabelled and "A"), respectively. 

The convergence is more complete at some points than at others. The 
non-convergance at z = 2 and beyond seems to be small -amplitude numerical 
"noise," but the over-all solution is sufficiently well defined for the 
present purposes. The difficulty of the convergence at and beyond z = 2 
may be associated in part with insufficiencies in the numbers of grid 
points, numbers of trajectories, and individual trajectory accuracy; 
impiovements in these parameters requires more computer time. At points 
along the axis more trajectories were used than at points off the axis, 
so that che acctiracy is relatively high along the axis. In contrast to 
off-axis points, the convergence is clearly excellent at the axial points 
beyond z = 1. The accuracy of n^. at off-axis points far downstream where 
({• is small is estimated to be about 1C percent. 

Tables 4-1 and 4-2 also give the dimensionless ion and electron cur- 
rent density (j. and j^) at the center of the wake side of the disk, (j^ 
is seen to have the value exp(-4).) The electron density (n^) profiles 
are given in Fig. 4-1 (middle profiles) and in Table 4-1 (in pa^'enthe'^cs).# 
n is assumed to be unity (ambient value) at the downstream boundary v: = 6 
where <{i is assumed to be zero. Figure 4-1 and Table 4-1 indicate that 
quasi neutral! tv (n = n^) is roughly valid within the accuracy of the cal- 
culation at t z = 4 and beyond. 

The features of the wake structure are as follcys: There is no prom- 
inent struv'ture such as large-amplitude ion or electron density bumps. 

The apparent structure in the ion profiles at z = 2 and beyond suggests 
the possibility of small -amplitude ion structure near the axis, beyond 
z = 2. The apparent ion structure near the axis at z - 5 is clearly not 
associated with the local potentidl profile (since f is essentially zero 
in this region); hence this structure must ba associated with upstream 
perturbations, i.e., the deflection of ion trajectories passing near ^he 
edge of the disk. It is also evident that the downstream boundary at z = 6 
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has been chosen sufficiently for downstream. This evidence is based on 
tte smoothness with which the potential has already fallen off to negli- 
gible values at z » 5. 

The reality of the small -aa^)1itude structures (as opposed to itera- 
tive noise) can be verified by more accurate calculations, with changes 
in the numerical parameters such as mmbers of grid points and trajec- 
tories. Persistence of the structure despite changes in the mmierical 
parameters may be taken as an indication of its reality. In spite of this 
uncertainty, the important gross features clearly indicated by the profiles 
are that (a) the wake becomes filled in by :ilectrons and ions sanewhere 
between 2 and 3 radii downstream, i.e., less than the Mach number of radii, 
and (b) the wake disturbance is essentially confined within a region 
extending to about z = 4 downstream, and outward to about r = 1.5 in the 
transverse dimension. 
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4,2 Case 4 q = -4, M * 4, and = 1/100 


For this case (Fig. 4-2 and Tables 4-4 to 4-6), the parameter values 
differ from those of the preceding problem only in x^, »diich is very small so 
that the problem applies to a large body, namely 100 Debye lengths in 
radius. This size of moving body is larger than has been previously 
treated by trajectory-following, i.e., realistic, calculations. (In the 
large-body calculations of Kiel et (1968), the particle trajectories 
were not treated realistically.) The results show what may be expected 
for th-' wake structure of large bodies in general. This case requires 
more effort (ccopjter time and judicious selection of numerical param- 
eters) than that of a smaller body. The solutions shown, therefore, are 
intended to be illustrative rather than accurate. 

Six iterations, or Poisson-Vlasov cycles, were computed in which 
successive iterates were used without mixing, starting with the neutral 
ion density as an initial guess. The nominal nunri>er of trajectories, 5i2, 
was used at all grid points. 

The profiles of n^, n^, and in Fig. 4-2 are constructed in the saro 
way and at the same grid points as in Fig. 4-1. The wake is essentially 
"empty" of both ions and electrons between z = 0 and z = 1, and begins to 
fill up between z = 2 and z = 3. In this way, the wake is qualitatively 
similar to that in Fig. 4-1. 

Again, two sets of ion-density profiles are shown on the left side 
of Fig. 4-2, the unlabelled profiles for the 6th order (6th iteration. 

Table 4-4), and the profiles labelled "A" for the 5th order (Table 4-5). 

The 6th-order potentials are given in Table 4-6, and the 6th-order n^- 
values are given, in parentheses, in Table 4-4. Comparison of the n^- 
values in Table 4-4 with the 5th-order n^-values (labelled "A") in Table 
4-5 indicate that the quasi neutrality assumption is valid everywhere out- 
side a cone-shaped region near the wake surface; the cone height along the 
axis is between one and two radii. This is in accord with expectation for 
a large body. Near the wake surface, however, quasineutrality is violated 
because the effective Debye length is large. The similarity of the n^- 
profiles (labelled "A") and the ng profiles in Fig. 4-2 is a consequence 
of near-quasineutrality. 
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Comparison of the 5th and 6th order n.-proflles (labelled "A* and 
unlabelled) in Fig. 4-2 show that the solution is reasonably converged 
for 2-1 and below, but that there is incomplete convergence at z = 2 
and beyond. The incomplete convergence and apparent structure at z = 2 
and beyond may be artifactual due to insufficient numerical accuracy. 

(No attempt was made to achieve high accuracy since this was regarded as 
a trial run.) The structure and lack of convergence are seen to extend 
past z = 5, so that the downstream boundary should be placed further than 
at 2 = 6. 

Coo|>ared with the previous case (and despite any inaccuracies), one 
may infer additional physical conclusions indicated by Fig. 4-2, nattily, 
(a) the suggestion of a core of high (approximately ambient) density of 
ions and electrons on the axis, and (b) the occurrence of a potential 
well in the near wake, defined as a region with ^-values below -4. The 
shading in the two lowest ♦-profiles of Fig. 4-2 denote cross-sections of 
this well. The wake-surface current densities (Table 4-4) are less than 
in the previous case; the electron current density is less than exp(-4), 
as would be expected in the presence of a potential well. 

The region of wake disturbance is not as well defined as in the pre- 
vious case of the smaller body, but it probably extends more than 6 radii 
downstream, and between 2 and 3 radii in fie transverse direction. 


- 27 - 



4.3 Case * -4, M = 8, and Xp = 1/100 


This case (Fig. 4-3 and Tables 4-7 to 4-9) is another large-body 
case similar to the previous large-body case except that the Mach ntxnber 
is inci'eased from M = 4 to H = 8. Ten iterations were computed in which 
successive iterates were used without mixing, starting with uniform 
ambient ion density. (The latter starting condition was inadvertently 
different from that of the M = 4 calculation which was started with the 
neutral ion density, but this difference should become unimportant after 
many iterations.) Similar statements may be made about the incomplete- 
ness of the convergence as in the M = 4 case. The 9th and lOth-order ion 
densities are labelled "A“ (Table 4-8) and unlabelled (Table 4-7), respec 
tively. On coR4)aring these, the convergence seems fairly good at z = 0.5 
and z = 1. Again, the disturbance extends beyond z = 5, so that the down 
stream boundary should be moved further than z = 6. 

Despite inaccuracies, the consistency is such that physical conclu- 
sions may be drawn as follows. In this case the wake is seen to remain 
anpty funther downstream than in the M = 4 case. In addition, the sugges 
tion is much stronger that there is a central core of ambient density for 
both ions and electrons along the axis. Moreover, the potential well is 
wider and longer than in the M = 4 case, although the depth is about the 
same. The disk-wake-surface electron current density (Tcble 4-7) is 
slightly less than the M = 4 value, and is again less than exp(-4). 

The conical region behind the disk where quasi neutrality breaks 
down is now longer than in the M = 4 case, extending to between z = 4 
and z = 5 along the axis. 

The region of wake disturbance is probably longer than 6 radii 
downstream, as in the M = 4 case, but may not extend beyond about 2 
radii in the transverse direction. 
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TABLE 4-1 


NORMALIZED ION (AND ELECTRON) DENSITY IN WAKE 
-4, M - 4, Xq ■ 1/5) 


z 

r=0 

r».l 

r«.3 

r«.6 

r«.8 

r-1.0 

r-1.2 

r-1.5 

r»2.0 

6 

1.03 

1.03 

1.09 

0.95 

0.95 

0.95 

1.09 

1.09 

1.09 


(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

5 

1.03 

0.79 

1.02 

0.95 

1.01 

1.02 

1.08 

1.06 

1.06 


(0.98) 

(0.98) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

4 

1.01 

0.89 

1.00 

0.99 

0.96 

0.97 

0.95 

0.95 

0.95 


(0.95) 

(0.95) 

(1.01) 

(1.02) 

(1.03) 

(1.03) 

(1.03) 

(1.03) 

(1.03) 

3 

0.96 

0.82 

1.01 

1.06 

1.07 

1.09 

1.09 

1.07 

1.07 


(0.91) 

(0.91) 

(0.97) 

(0.99) 

(0.99) 

(0.99) 

(0.99) 

(0.99) 

(0.99) 

2 

0.39 

0.46 

0.53 

0.51 

0.64 

0.73 

0.79 

0.96 

0.95 


(0.58) 

(0.58) 

(0.61) 

(0.63) 

(0.69) 

(0.75) 

(0.83) 

(0.96) 

(1.00) 

1 

0.016 

0.018 

0.016 

0.15 

0.49 

0.68 

0.79 

0.88 

0.98 


(0.19) 

(0.19) 

(0.21) 

(0.31) 

(0.45) 

(0.61) 

(0.74) 

(0.87) 

(0.96) 

0.5 

0.00034 

0.00037 

0.00047 

0.0035 

0.20 

0.77 

0.91 

0.88 

0.97 


(0.084) 

(0.085) 

(0.094) 

(0.14) 

(0.23) 

(0.42) 

(0.65) 

(0.84) 

(0.95) 

0.2 

2.3x10'^ 

2.2x10"® 

2.6x10"^ 

2.0x10"® 

0.00039 

0.82 

1.00 

0.93 

0.97 


(0.036) 

(0.036) 

(0.038) 

(0.049) 

(0.071) 

(0.15) 

(0.46) 

(0.81) 

(0.94) 

r, z 

in units of disk radius 




IT - 12, 9/8/75, . 

a ■ .5 

ji = 

9.6x10"®, 

jg = 0.0183 





(j) ■ 0 

at 2 „ 




TABLE 4-2 


NORMALIZED ION DENSITY IN WAKE (A) 
(♦q = -4, M ■ 4, Ap » 1/5) 


z 

r=0 

r=.l 

r=.3 

r*.6 

r=.8 

r=1.0 

r-1.2 

r«l .5 

r*2.0 

6 

1.04 

1.04 

0.95 

0.95 

0.95 

0.95 

0.95 

0.95 

0.95 

5 

1.04 

0.91 

1.00 

1.00 

1.00 

0.99 

0.97 

0.97 

0.97 

4 

1.00 

0.77 

1.05 

1.05 

1.08 

1.06 

1.09 

1.07 

1.07 

3 

0.96 

0.59 

0.99 

0.99 

0.98 

0.97 

0.95 

0.96 

0.94 

2 

0.40 

0.45 

0.71 

0.54 

0.60 

0.72 

0.79 

0.98 

1.03 

1 

0.017 

0.018 

0.016 

0.15 

0.50 

0.67 

0.79 

0.88 

0.98 

0.5 

0.00035 

0.00037 

0.00048 

0.0035 

0.19 

0.77 

0.91 

0.90 

0.98 

0.2 

0.000022 

0.000022 

2.6x10"^ 

0.000020 

0.00039 

0.82 

1.00 

0.93 

0.97 


9.4x10'®, 

jg * 0.0183 





IT = 

n. 9/8/75, 

a * .5 


♦ = 0 at 
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TABLE 4-3 


NORMALIZED POTENTIAL IN WAKE 
(♦q- -4. M » 4» Xjj » 1/5) 


z 

r =0 

r*.l 

r *,3 

ra .6 

r *.8 

r » 1.0 

r*l ,2 

r » 1.5 

r = 2.0 

6 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

5 

- 0.016 

- 0.019 

- 0.00046 

0.0016 

0.0016 

7 . 5 x 10 '® 

- 0.0037 

- 0.0041 

- 0.0027 

4 

- 0.046 

- 0.048 

0.0075 

0.022 

0.030 

0.029 

0.033 

0.030 

0.026 

3 

- 0.096 

- 0.099 

- 0.027 

- 0.011 

- 0.0073 

- 0.0058 

- 0.0064 

- 0.0053 

- 0.0083 

2 

- 0.54 

- 0.54 

- 0.49 

- 0.47 

- 0.37 

- 0.28 

- 0,18 

- 0.037 

0.0013 

1 

- 1.66 

- 1.64 

- 1.54 

- 1.17 

- 0.79 

- 0.50 

- 0.30 

- 0.14 

- 0.045 

0.5 

- 2.48 ' 

- 2.46 

- 2.36 

- 1.98 

- 1.49 

- 0.86 

- 0.42 

- 0.17 

- 0.055 

0.2 

- 3.33 

- 3.32 

- 3.26 

- 3.03 

- 2.65 

- 1.87 

- 0.78 

- 0.21 

- 0.059 


C ::d 



TABLE 4-4 


NORMALIZED ION (AND ELECTRON) DENSITY IN WAKE 
(4b= -4. M = 4, Xjj = 1/100) 


2 

r=0 

r=.l 

r=.3 

r=.6 

r=.8 

r=1.0 

r»1.2 

r«l .5 

r«2,0 

6 

1.08 

(1.00) 

0.87 

(1.00) 

0.87 

(1.00) 

0.80 

(1.00) 

0.88 

(1.00) 

0.91 

(1.00) 

0.98 

(1.00) 

0.91 

(1.00) 

0.92 

(1.00) 

5 

1.09 

(0.91) 

0.79 

(0.52) 

0.75 

(0.85) 

0.78 

(0.75) 

0.95 

(0.75) 

0.91 

(0.84) 

0.97 

(0.89) 

0.90 

(1.07) 

0.91 

(1.06) 

4 

1.00 

(0.81) 

0.28 

(0.76) 

0.58 

(0.63) 

0.76 

(0.65) 

0.95 

(0.78) 

0.85 

(0.79) 

0.81 

(0.81) 

0.89 

(0.91) 

0.99 

(0.94) 

3 

0.76 

(0.52) 

1.12 

(0.50) 

0.65 

(0.66) 

0.72 

(0.62) 

0.63 

(0.60) 

0.75 

(0.72) 

0.81 

(0.80) 

0.88 

(0.85) 

0.95 

(0.94) 

2 

0.96 

(0.57) 

0.24 

(0.65) 

0.44 

(0.63) 

0.36 

(0.42) 

0.36 

(0.43) 

0.67 

(0.67) 

0.76 

(0.76) 

0.85 

(0.87) 

0.95 

(0.95) 

1 

0.14 

(0.044) 

0.12 

(0.13) 

0.039 

(0.049) 

0.11 

(0.16) 

0.30 

(0.27) 

0.59 

(0.64) 

0.75 

(0.74) 

0.86 

(0.86) 

0.96 

(0.97) 

0.5 

0.00020 

(0.0047) 

0.00093 

(0.0050) 

0.0037 

(0.0063) 

0.019 

(0.024) 

0.13 

(0.19) 

0.57 

(0.55) 

0.79 

(0.73) 

0.84 

(0.89) 

0.94 

(0.97) 

0.2 

0.0095 

(0.0038) 

2.7x10'^ 

(0.0038) 

1.3x10"^ 

(0.0041) 

2.8x10"^ 

(0.0059) 

0.00054 

(0.013) 

0.56 

(0.56) 

0.80 

(0.75) 

0.83 

(0.91) 

0.94 

(0.98) 


2.4x10'^, 

= 0.0043 
e 
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TABLE 4-5 


M0W-1ALIZED ION DENSITY IN 
(<t>0= -4. M = 4, Xp » 1/ 


z 

r=0 

r=.l 

r=.3 

r=.6 

r».8 

6 

0.94 

0.84 

0.84 

0.77 

0.86 

5 

0.94 

0.52 

0.86 

0.75 

0.75 

4 

O.Bl 

0.76 

0.63 

0.65 

0.78 

3 

0.52 

0.00 

0.66 

0.62 

0.60 

2 

0.56 

0.65 

0.63 

0.42 

0.43 

1 

0.00063 

0.14 

0.046 

0.16 

0.27 

0.5 

0.0017 

0.(025 

0.0033 

0.021 

0.19 

0.2 

0.000035 

5.4x10"^ 

8.9x10'® 

0.000033 

0.00050 


= 1.1x10"®, jg = 0.0044 


(A) 


r=1.0 

r»1.2 

r*l .5 

r=2.0 

0.96 

0.94 

1.09 

1.09 

0.84 

0.89 

1.07 

1.06 

0.79 

0.81 

0.91 

0.94 

0.72 

0.80 

0.85 

0.94 

0.67 

0.76 

0.87 

0.95 

0.64 

0.74 

0.86 

0.97 

0.56 

0.73 

0.89 

0.97 

0.57 

0.75 

0.91 

0.98 
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TABLE 4-b 


NORMALIZED PO'^ENTIAL IN WAKE 
(4,^= .4, M = 4. Xp = 1/100) 


z 

r=0 

r=.l 

r=.3 

r=.6 

r».8 

r=1.0 

r*1.2 

r=1 .5 

r»2.0 

6 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

5 

-0.090 

-0.65 

-0.16 

-0.29 

-0.29 

-0.17 

-0.12 

0.065 

0.057 

4 

-0.21 

-0.28 

-0.46 

-0.44 

-0.24 

-0.24 

-0.21 

-0.095 

-0.057 

3 

-0.65 

-0.68 

-0.42 

-0.48 

•0.51 

-0.33 

-0.22 

-0.17 

-0.058 

2 

-0.56 

-0.44 

-0.46 

-0.87 

-0.83 

-0.40 

-0.28 

-0.14 

-0.047 

1 

-3.12 

-2.03 

-3.01 

-1.86 

-1.30 

-0.45 

-0.30 

-0.15 

-0.031 

0.5 

-5.35 

-5.30 

-5.07 

-3.72 

-1.67 

-0.59 

-0.32 

-0.12 

-0.032 

0.2 

-5.57 

-5.56 

-5.49 

-5.14 

-4.32 

-0.58 

-0.29 

-0.095 

-0.022 


TABLE 4-7 


NORMALIZED ION (AND ELECTRON) DENSITY IN WAKE 
-4, M = 8, Xp = 1/100) 


z 

r=0 

r=.l 

r=.3 

r=.6 

r=.8 

r=1.0 

r=1.2 

r=l .5 

r»2.0 

6 

1.11 

1.09 

0.71 

0.80 

1.01 

1.32 

1.10 

1.10 

1.11 


(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(1.00) 

(l.CO) 

5 

1.11 

1.24 

0.045 

0.80 

0.63 

0.91 

1.12 

1.11 

1.12 


(1.04) 

(0.18) 

(0.038) 

(0.44) 

(1.30) 

(0.90) 

(1.10) 

(1.11) 

(1.12) 

4 

1,04 

0.0019 

0.24 

0.63 

0.64 

0.96 

1.04 

1.14 

1.15 


(1.21) 

(0.080 

(0.057) 

(0.38) 

(0.54) 

(0.87) 

(1.20) 

(1.13) 

(1.14) 

3 

1.27 

0.015 

0.095 

0.0065 

0.71 

1.08 

0.99 

1.24 

1.27 


(0.93) 

(0.055) 

(0.31) 

(0.023) 

(0.37) 

(0.95) 

(0.99) 

(1.13) 

(1.10) 

2 

1.04 

0.00064 

0.00081 

0.0070 

0.46 

0.85 

0.98 

1.12 

■ 1.06 


(1.12) 

(0.022) 

(0.11) 

(0.061 ) 

(0.47) 

(0.85) 

(1.05) 

(1.26) 

(1.21) 

1 

1.27 

1.5x10'^ 

9.7x10"^^ 

0.00019 

0.0014 

0.87 

1.01 

1.28 

1.32 


(0.85) 

(0.0071) 

(0.0025) 

(0.0039) 

(0.014) 

(0.79) 

(0.94) 

(1.07) 

)1.08) 

0.5 

5.0x10"^^ 

4.4x10'^^ 

3.0x10"^ 

4.6x10"^^ 

8.6x10'^ 

0.90 

1.01 

1.12 

1.08 


(0.0018) 

(0.0017) 

(0.0018) 

(0.0034) 

(0.013) 

(0.82) 

(0.95) 

(1.11) 

(1.21) 

0.2 

6.9x10"^ 

1.6x10"^ 

9.4x10'^® 

1.3x10"^° 

3.6x10"® 

0.87 

1.16 

1.18 

1.18 


(0.0032) 

(0.0032) 

(0.0036) 

(0.00/7) 

(0.20) 

(0.86) 

(0.96) 

(0.93) 

(1.30) 


4.2x10'^°, 

= 0.0037 
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TABLE 4-8 


NORMALIZED ION DENSITY IN WAKE (A) 
(♦q- -4, M - 8. Xp - 1/100) 


z 

r-*0 

r-.l 

r-.3 

r",6 

r-.8 

r-1.0 

r-1.2 

r-1.5 

r-2.0 

6 

1.11 

1.07 

0.33 

0.79 

1.33 

1.33 

1.10 

1.10 

1.11 

5 

i.n 

0.18 

0.032 

0.44 

1.31 

0.90 

1.10 

1.11 

1.12 

4 

1.31 

0.075 

0.052 

0.38 

0.54 

0.87 

1.20 

1.13 

1.14 

3 

1.04 

0.039 

0.31 

0.011 

0.87 

0.95 

0.99 

1.13 

1.10 

2 

1.27 

0.0038 

0.11 

0.056 

0.47 

0.86 

1.05 

1.26 

1.21 

1 

1.04 

0.030017 

O.OOOU37 

0.00095 

0.0044 

0.80 

0.94 

1.07 

1.08 

0.5 

2.1x10"^^ 

4.6x10’^® 

3.0x10’® 

3.8x10’® 

0.0014 

0.83 

0.95 

1.11 

1.21 

0.2 

6.9x10"^ 

6.2x10“^ 

4.3x10’^ 

8.2x10’^ 

0.21 

0.87 

0.96 

0.93 

1.31 


7.4x10"^°, 

jg - 0.0029 
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TABLE 4-9 


MORMALIZEO POTENTIAL IN WAKE 
(4b“ -4, M ■ 8. jj » 1/100) 


z 

r=0 

r=.l 

r-.3 

r«.6 

r».8 

r-1.0 

r-1.2 

r-1.5 

r«2,0 

6 

0.0 

0.0 

0.0 

0.0 

0.0 

O.C 

0.0 

0.0 

0.0 

5 

0.039 

-1.71 

-3.26 

-0.82 

0.26 

-0.10 

0.094 

0.11 

0.12 

4 

0.19 

-2.52 

-2.87 

-0.96 

-0.62 

-0.14 

0.19 

0.13 

0.13 

3 

-0.073 

-2.90 

-1.18 

-3.76 

-0.14 

-0.047 

-0.015 

0.12 

0.094 

2 

0.11 

-3.80 

-2.24 

-2.79 

-0.75 

-0,16 

0.051 

0.23 

0.19 

1 

-0.16 

-4.95 

-5.98 

-5.55 

-4.30 

-0.23 

-0.059 

0.070 

0.081 

0.5 

-6.33 

-6.38 

-6.31 

-5.69 

-4.38 

-0.19 

-0.053 

0.11 

0.19 

0.2 

-5.75 

-5.74 

-5.63 

-4,86 

-1.60 

-0.15 

-0.044 

-0.077 

0,27 


APPENDIX A 


THE VLASOV PROBLEM: DENSITY CALCULATION 

For the purpose of evaluating density and current-density i;itegrals, 
it is convenient to transform to energy and angle variables in velocity 
space. Since v»e will be interested primarily in Maxwellian energy distri- 
butions (with drift), we adopt the following units in terms of which dimen- 
sionless variables may be defined: 

kT = unit of energy, where T is the temperature of the 
Maxwellian distribution 

v^lHTra = unit of velocity, namely, the most probable velocity 

n^ = unit of particle density, nan^ly, the unperturbed 
density 

The energy and angle variables are: 

E = energy in muHiples kT 

a = polar angle with respect to z-axis (Fig. A-1) 

B = azimuthal angle with respect to the plane containing 
the z-axis and the point r (Fig. A-1) 

z-axis = axis of sy.naietry of body as well as direction of 
plasma flow 

The definitions of the angles of a and 6, which define the orientation of 
the velocity-vector v, are illustrated in Fig. A-1. The potential energy 
t will also be a multiple of kT. 
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The integral for the particle density follows from appropriate trans- 
formation of the volume element in velocity space of Chapter 2, namely. 


-111 


f V dv sin a da d8 


= (const) 


111 '■ 


/E - ^ dE sin o da d8 


(A-1) 


where there is a step-function s (a "cut-off" factor) ihicli is defined as 
unity if the orbit connects with infinity (an "escaping" orbit), and as 
zero otherwise. The step-function thus automatical ly takes care of the 
restriction to the domain of escaping orbits. Hovjever, whether 6 is unity 
or zero is decided only after the orbit has been followed backward in tiim? 
sufficiently far (by performing an orbit calculation in the given potential 
distribution). 

For the Maxwellian distribution with drift, v;e have 


n = 


2n 


3/2 


dE 


sin a drt 


2it 

^ -U(E,a,6) 


Max (0,^) 


SdS 


(A-2) 


where 


U = E + + 2M^ cosa (A-3) 

CD 

with H denoting the Mach number (plasma-flow velocity divided by /21t/m) , 
and denoting the value of the polar angle a of the velocity-vector at 

infinity. 

The limits on the integral correspond to the full ranges of the vari- 
ables. The lower limit, "Max(0,<ji)," on the E-integral is defined to be 
zero if (}> < 0 (attractive potential), and to be if t> > 0 (repulsive 
potential). For = 0 and 6 = 1, is equal to a and n becomes unity (the 
unperturbed value). 
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The current density integral may similarly be written 


j * 




dv cos a sin a da dg 


* (const) 



(E - ♦)dE cos a sin o da d3 


(A-3) 


where the normal component of velocity is identified with the 
2 «c(xnponent of velocity, if the surface of interest on which the current 
is to be calculated ^s perpendicular to the z-axis. The constant in 
front will now be adjusted. Specializing Eq. (A-5) to the Maxwellian distri- 
buticn (with drift), we may write 


Jo » J.. 


r 

(E - ♦)JE cos a sin a da 


Max(0,^) 


■2x 


j-U(E,a,e) 


The current density is expressed as the ratio j/j^ where denotes the 
"random" current density (namely, n^.^72iTm) which would be collected in 
the absence of plasma flow and if there were no electric fields, i.e., 
such that ^ = 0 and 6=1 over the ranges of integration. In the remainder 
of the report it will be convenient to take j to mean j/j^ (j^ 3 1). 

The integrals in Eqs. (A-2) and (A-4) represent infinite numbers of 
orbits (continuous distributions in a, ^nd g). In approximating the 
integrals by quadratures, we replace the infinite sums by finite numbers 
of terms, each corresponding to an orbit. Thus, transfurming the ranges 
of integration into intervals between -1 and +1 in preparation for the 
Gaussian quadrature, we may write 

E(c) = Max(0,^) (A-5) 
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a 


a 


3 


cos"^ a 


for density 
sin"^^ for current density 

fo^b) 


(A-6) 


(A-7) 


Then the transfomed density and current density integrals become 


and 


,1 


n = A 


1 


3 = 0 


-U(c,a,b) 


/ECc) - ^ * 6 • 


n n 


dc da db 

(1 - c)2 


(A-8) 


rl 


.1 


g-U(c,a,b)|-^^^j 


- ♦] * « 


-1 -1 -1 


dc da db 
(1 - c? 


(A-9) 


In a case where a potential barrier exists, "Max (0,<^)" in Eq. (A-5) 
should be replaced by the barrier height. 


We now have the integrals in a form suitable for Gaussian quadratures, 
where the new variables (c, a, b) all lie in the range -1 to +1. For flex- 
ibility, we now divide the c-range into sub- intervals, and apply a 
Gaussian quadrature of order 2 to each of these sub-intervals. Similarly, 
we divide the a-range into sub-intervals, and the b-range into sub- 
intervals, with Gaussian quadratures of order 2 in each sub-interval. Then 
both Eqs. (A-8) and (A-9) may be put in the form 


rl rl rl 


I = 


T(E) • 6(E, a, e) • dc da db 


1 ll 


(A-10) 


which may be approximated by the sum: 
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(A-11) 


(A-12) 


with 

E' H E(c’), E” H E(c") 

a’=o(a'), a” = a{a”) (A-13) 

e'=e(b'). 3" = 0 (b") 

and with E(c), a(a) and g(b) defined by Eqs. (A-5) through (A-7). 

The special values c‘, c", a', a", b', and b" are defined by formulas 

based on the abscissas for the Gaussian (Order-2) quadrature 

applied to the multiple sub-interva's, namely. 
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c* = 


1 


- + 2K - 1 - M 
/I ® ® 


c" = — 


+ — + 2K - 1 - 
✓J ® ® 


a' = 


1 


- — + 2K - 1 
/I * 


- M. 


1_ 

M. 


+ — + 2K - 1 
/3 ^ 


- M. 


b‘ = i- 




- — + 2K^ - 1 
/3 '' 


- M. 


= it 


+ — + 2K. - 

/3 ^ 


1 - M. 


(A-14) 


Again, the function 6 is the unit step-function which is unity for escap 
ing orbits. Note that the Gaussian coefficients are (conveniently) all 
unity for this quadrature scheme. 

In the absence of plasma flow (M = 0), or for electrons, one can 
consider I to be a sum of monoenergetic contributions, which becomes evi 
dent by rearranging the sum, namely: 



[T(E') • F(E') + T(E") • F(E")] 


(A-15) 


where F(E) is a new quantity defined by: 


F(E) H 


i Mb 

kZI 

^=1 V ' 


[6(E. a'. S') + 6(E, a", S")] (A-16) 
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Note that, since 5 is either unity or zero according as the orbit escapes 
or is absorbed, the sum of Eq. (A-16) is the quotient of two integers, 
namely, the number of orbits escaping divided by the total number of 
orbits for the given energy E. Thus, F is a fraction between 

zero and unity, and becomes unity if all orbits escape. Note also that 
information regarding the energy distribution resides in T. Thus, a non- 
Maxwellian distribution may be treated by suitably modifying T, In par- 
ticular, for a monoenergetic distribution we simply set E equal to unity 
and replace S, Eq. (A-15), by the single term: 


S 


mono 



• F(l) for 


{ current density 
density 


(A-17) 


where F(l) is given by Eq. (A-16) with E = 1 and represents the fraction 
of orbits which escape. The dimensionless potential <j> is now a multiple 
of the singular energy of the particles. 


The equations derived here are suitable for a computer program and 
have been incorporated into the program used for the results discussed 
in Chapter 4. 


The method of computation of orbits involves integration of the 
equations of motion, with the forces given by the components of the gradi- 
ent of potential. These components are obtained by interpolation between 
values of potential defined at the points of a grid in space as described 
in the next appendix. The criterion for "escape" of an orbit (i.e., evalu- 
ation of 6) depends on the geometry of the problem and of the grid. The 
equations of motion are integrated step-by-step until the orbit either 
passes out of the outer boundary of the grid ("escapes" so that 6 = 1) or 
returns to one of the metal surfaces (is "absorbed" so that 5 = 0). The 
orbit computation time-step is not of physical importance in these time- 
independent problems where only the shape of the orbit matters. The time- 
step is kept as large as possible consistent with maintaining the energy 
loss or gain within desired limits. The method of integrating the equations 
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of motion, the interpolation method to find the forces, and the control 
of step size, are discussed in Appendix C. 

An important consideration is the accuracy of the quadrature-sum. 
Naturally, the accuracy is related to the number of terms used, that is, 
the number of orbits where each term corresponds to a unique orbit. The 
total number of orbits involved in Eq. (A-11) or Eq. (A-15) is given by 
8M M H. . In a test of the energy quadrature alone, and with M = 0, the 
unperturbed value of density (unity) vias computed for values of = 1, 2, 

4, 8, 16, and 32. The corresponcing numerical errors were -6», -7%, 

+1.5%, -0.05%, +0.013%, and +0.0D3?.. This test was independent of geom- 
etry (the o and B integrations were numerically exact). Thus, = 4 
(8 values of E) is taken to represent sufficient accuracy (within a few 
percent) for the purposes of computing density for a Maxwellian distri- 
bution without drift (or, for electrons). 

A device for improving the accuracy of the quadratures at large Mach 
numbers, without increasing the total number of orbits and therefore com- 
puter time, is to suitably weight the integrand in the domains of importance. 
Thus, one modification is to multiply the term "(I + c)/(l - c)" in Eq. 

(A-5) by M when M exceeds a suitable value, say unity. This emphasizes the 
contributions of orbits having E in the vicinity of M . Another modifica- 
tion is to replace "a" in Eq. (A-6) by the function "-1 + 2((1 + a)/2)^“. 

This emphasizes the contributions of orbits having a in the vicinity of ir, 
that is, directed upstream. When these modifications are used, the quad- 
rature sums must also be multiplied by additional corresponding factors, 
namely, and "M((l + a)/2)^^'^", respectively. 
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APPENDIX B 


THE POISSON PROBLEM: POISSON DIFFERENCE EQUATIONS 

In the problems treated in this report the electrostatic field is 
axially symmetric and is defined on a mesh of spatial grid points, such 
that at any point (including grid points) the potential and electric 
field are obtained by interpolation. 

Assume that the space charge density is known at the grid points. 
Consider a group of interior grid points, forming a portion of the over- 
all grid as shown in Fig. B-1. In this figure, the vertical and hori- 

zontal directions are the z and r directions, respectively, where z and 
r denote the cylindrical axial and cylindrical radial coordinates, 
respectively. Three horizontal grid lines, of constant z-values 
z^, and and three vertical grid lines, of constant r-values 

rj, and r^^^, are shown in the figure. (Note that the index (i) of z 
increases as z decreases.) ihe set of grid lines intersect at 9 grid 
points, or nodes, as shown. Each point may be considered to be asso- 
ciated with a volume of space, and to have a group of four neighboring 
points which "interact'' with it. Thus, consider the central point of 
the group, labelled C in the figure. Associated with this point is a 
volume of revolution (a torus) whose r~oss-section is rectangular and 
is shown by the rectangular shaded area surrounding Point C. The shaded 
area is defined by connecting the mid-points of the surrounding mesh 
rectangles. Let x denote the volume of the torus, and let the neighbor- 
ing points (above, below, to the right of, and to the left of C) be 
labelled N, S, E and W (north, south, east and west, respectively). 

Let the Poisson equation be written in dimensionless form as 


= 




ni)/x jj 


(B-1) 


where n^, n^, 41 and p denote the dimensionless electron density, ion 
density, Debye length, electrostatic potential and space-charge density, 
respectively; and all lengths are in units of the body radius. Now inte- 
grate Eq. (B-1) over the volume x of the torus associated with point C: 
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where Pq is known at the grid point C. The right-hand side has been 
approximated as shown since t is small in principle, and Pq is the value 
of p at Point C. By the divergence theorem, the left-hand side becomes 



v.'hera z denotes the surface of the torus; 94>/3n is the component of 
in the outward normal direction at the surface; A^, A^, and A^^ denote 
the areas of the north, south, east, and west surfaces, respectively; and 
the quantities ( 3 (}>/ 3 n)j^ SEW values of d<p/^n taken to be constant 

on the corresponding surfaces. 

( 3 (j)/Dn)^ SEW approximated by difference quotients, namely. 




34>1 * ^*W " 


(B-5) 

(B-6) 

(B-7) 
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where # denotes the potential at Point C and denote tt^ 

neighboring potentials. The areas Aj^, A^, A^, and Ay arc given by 


S = J KN*, ♦ rj)* - (rj ♦ rj_,)2] 

(B-8) 

II 

(B-9) 

*E ' 7 ■ ^i»l- 

(B-10) 

*« ' f f'^j * '■j-lX^i-1 - 

(B-11) 

and the volune t is given by 


Ay 

T = 2“ ^^i-1 ' ^i+1^ 

(B-12) 

Thus, equating Eq. (8-3) with Eq. (B-2), and substituting 
we obtain the difference equation in the form 

the foregoing. 

S *N * ^sH * ^E^E * ~ ^ ~ Pc^ 

(B-13) 

where 



(B-14) 

and 



(B-15) 

*5 

(B-i6) 
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(B-17) 






(B-18) 


This shows how to form the difference equations used for the Poisson prob- 
lems of this report. Equation (d-14) holds only for an "interior” point 
of the grid, that is, a >x)int surrounded by neighbors on all four sides. 

If Point C 'S not an interior point, one or more of the terms on the 
left-hand side of Eq. (B-13) may vanish. To see how this happens, con- 
sider an exanyile griH shown in Fig. B-2 where there are 16 grid points, 
but only 13 of these correspond to unknown potentials to be solved for.* 

The grid points where the potentials are unknown are numbered and indi- 
cated by circles. The three solid circles labelled by the letters a, b, 
and c denote electrode sur<"ace ooints where the potentials are known. 

Points No. 1, 2, 4, 5, 9, 10, and 12 are special points all of vvhich have 
different properties, and are indicated by small crosses wihin the cir- 
cles. Among these points, the only interior point is No. 10. 

Consider Point No. 10, which has a known potential, namely = <}i^, 
as its southern neighbor. The equation for this point is given by a modi- 
fication of Eq. (B-13) in that the term is now known and is trans- 
ferred to the right-hand side. For this point, Eq. (B-14) still holds. 

Consider Point No. 5, which is on the axis. Its equation is given 
by the modification of Eq. (B-13) in which Cy vanishes, and the remaining 
coefficients are evaluated with r. - r. . = 0. Equation (B-14) still holds, 
but with C^ net appearing. 

Consider Point No. 9, which is both on the axis and has a known neigh- 
bor (<)^ = 1 ^^), The modification of Eq. (B-13) includes both the modifica- 
tions for Points No. 5 and 10. Equation (B-14) still holds, but with C^ 
not appearing. 


♦This grid illustrates only the space behind the body; at points else- 
where around the body the formulas are similar. 
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North boundary/ of grid 



surface points 

FIG. B-2. SPECIAL GRID POINTS 
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East boundary of grid 


The modifications for Points No. 2* 4, and 12 depend on the fact 
that these points are on the outer boundary of the grid, where a “float- 
ing" boundary condition is used. Namely, the potentials are assumed to 
satisfy the linear laws 


11 = 11 = 


an 3Z 


A'4 


(B-19) 


and 


11 = 11 = 
an ar 


B’4 


(B-20) 


on the z-boundary (north, or south), and on the r-boundary (east), re- 
spectively. The formulas chosen for A* and B' depend on the assumed 
potential model. As a result of using Eq. (B-19) or (B-20) (or both for 
a comer point), corresponding terms do not appear on the left-hand side 
of Eq. (B-13), and the coefficient C is suitably modified so that Eq. 
(B-14) no longer holds. That is, for Points 1, 2, 4, and 12, C is given. 


respectively, by 

C = A‘Aj^ + €5 + 0 ^ (Point 1) (B-21) 
C = A'Aj^ ^ S ^ ^E ^ 2) (B-22) 
C = A'A|J + + B*A^ + Cy (Point 4) (B-23) 
C = C|^ + Cj + B'A^ + Cy (Point 12) (B-24) 


where it is understood that the coefficients in the above equations depend 
on the location of the point in question. (Aj^ and A^ denote modifications 
of A^j and A^; see examples given in Parker (1968).) 
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Once the coefficients of all of the equations (corresponding to the 
grid points where the potentials are unknown) are computed, the system of 
linear equations may be solved by iteration. Point-successive over- 
relaxation has been found to work well (Parker, 1968). 

The relaxation procedure may be illustrated by re-writing Eq. (B-13) 
in the form 

+ Cy#y j”* ” ^ (B-25) 

The superscripts n and n+1 denote the n-th and (n+l)-th iterates, respect- 
ively. At the beginning of each "sweep" through the equations, all are 
considered known. Then new values of ^ are obtained by putting on the 
right-hano side of Eq. (B-25) the most recently updated values of at the 
neighboring grid points. The superscripts “n, n+1" on the right-hand side 
of Eq. (B-25) ini^jly that one or more of the quantities <}>j, or 
might have been already updated by appearing on the left-hand side of a 
previous equation during the sweep. Thus far, the iteration sch«ne indi- 
cated by Eq. (B-25) is known as the "Gauss-Seidel" iteration. Convergence 
usually requires many iterations for problems involving many grid points. 

A more efficient procedure which will reduce the number of iterations 
is to modify the obtained from Eq. (B-25) before going on to the next 
equation in the sweep. The modification is given by the "relaxation" equa- 
tion 


^♦"^'Udified = “^*"^^^6auss- ^ 

Seidel 

where u is a relaxation parameter and the first term on the right-hand 
side of Eq. (B-26) involves as a factor the resulting from Eq. (B-25). 
The paranreter u may be less than unity ("under-relaxation") or greater than 
unity ("over-relaxation"). In practice it is found that the number of iter- 
ations is minimized dramatically when w is close to but less than 2. For 
example, simply increasing the value of w from 1.0 (Gauss-Seidel iteration) 
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to 1.9 (over-relaxation) typically reduces the numbers of iterations, in 
the problems of this report, from the order of thousands to the order of 
hundreds (for convergence to within a fractional change of one part in 
10® in all potentials). 

A single modification of the foregoing procedure "forces" the poten- 
tial to have arbitrary desired values at any of the grid points. This 
means that one may, for example, set the potential to zero at the outer 
grid points if for some reason this is felt to give a better representa- 
tion of "infinity" than the floating one. Or it allows one to arbitrarily 
specify the shape of the body surface, which may be of value for three- 
dimensional problems. The modification consists of re-defining (over- 
riding) the appropriate values of potential as soon as they are updated 
on the left-hand side of Eq. (B-25). This procedure is valid since the 
convergence of the iteration is not significantly affected. 

Another modification which is of value for large bodies is afforded 
when the potential distribution is such that is approximable by the 

Boltzmann factor exp(4). In this case, one may replace Eq. (B-25) by 
the equation 

+ (t/X^) exp(#"^M = + 03.^5 + 


+ n.T/xg (8-27) 

and then solve this nonlinear equation for with n. corsidered to be 

held fixed. 



APPtNDIX C 


COMPUTER PROGRAM 

This program Implements the calculational procedures described in 
Chapters 2 and 3, and in Appendices A and B. The main controlling pro- 
gram name is DWAKE (for "Disk Wake") or TDWAKE (for "Thick-Disk Wake," 
referring to modifications for a thick disk or short cylinder). There 
are two principal subprograms (FIELD and DENSTY), and eight subroutines 
serving these programs. The control hierarchy is listed as follows. 

DWAKE or TDWAKE (main) 

FiELO (Poisson problem) D ENSTY (Vlasov problem) 

LEFT INTERP 

MIDDLE TPJ\CK 

RIGHT 

PRINT 

SEIDEL 

ROOT 

The program operation is as follows. 

FIELD solves the Poisson problem for the potentials v;hen (optionally) 
either the ion or charge densities are given. DENSTY solves the Vlasov 
problem for the number and current densities vihen the potentials are given. 
The subroutines LEFT, MIDDLE, RIGHT, PRINT, and SEIDEL are called by FIELD, 
while INTERP and TRACK are called by DENSTY. FIELO and its subroutines con 
stitute one overlay, while DENSTY and its subroutines constitute the other 
overlay. SEIDEL in turn calls ROOT if the ion-dnsity option is selected. 

DWAKE reads the input data (described below) consisting of toe posi- 
tions of the grid lines, the values of the potential at points on the metal 
surfaces, the Debye number, the Poisson-Vlasov major-iteration mixing 
parameter a, the maximum number of majoi iterations, the array of input ion 
or charge densities, and various parameters affecting the options, orbit 
calculations, quadrature orders, single space-point coordinates, and single 
orbit initial conditions. 
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A Poisson-Vlasov iteration cycling loop is set up in which DWAKE 
calls FIELD and DENSTY, in each iteration cycle. DENSTY is called twice, 
once to compute the ion or charge densities and then to compute the cur- 
rent density. (In the appended listing the current density is computed 
only at the central wake point.) The cycling is continued for the desired 
maximum mrnibe of iterations, with the potentials or densities (normally, 
the potentials) being punched out on cards at every cycle for possible use 
in continuation. 

C.l Operation of FIELD and Its Subroutines 

FIELD sets up the coefficients of the potential -unknowns in the linear 
equations resulting from the differencing of the Poisson equation as des- 
cribed in Appendix B. Subroutine LEFT is called if the point in question 
is on the axis; subroutine RIGHT is called if the point is on the right- 
hand grid-boundary line; subroutine MIDDLE deals with an interior point. 

For each point corresponding to an unknown value of potential, FIELD obtains 
the values of the variables C, CN, CS, CE, and CW corresponding to the coef- 
ficients used in the equations of Appendix B. The values of these coeffi- 
cients, as well as the volume of the volume element (V) associated with the 
point in question, are printed out and saved (in the order in which they 
would appear as matrix elements) by subroutine PRINT. The right-hand sides 
of the Poisson difference equations are also set up u^ing the input ion or 
space-charge densities. When the setting-up has been completed, FIELD calls 
SEIDEL to obtain the solution, which is accomplished by point-successive 
over-relaxation. If the option is selected in which the nonlinear equation 
at the end of Appendix B is to be solved (appropriate for large bodies and 
where the Boltzmann factor can be used for the electrons), then SEIDEL calls 
ROOT within its over-relaxation loop. In ROOT, solution of each nonlinear 
equation is achieved by Newton's method for root evaluation. If the thick- 
disk option is selected, the body potential overrides the calculated poten- 
tials at grid points covered by the body. After SEIDEL has obtained the 
potentials in the form of a one-dimensional solution-vector X, the main pro- 
gram converts these to the two-dimensional potential arrays PHIN and PHIS 
used by DENSTY for the Vlasov problem. 
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C.2 Operation of OtNSTY and Its Subroutines 

OENSTY sets up an outermost loop in which one or more space points 
(or all grid points) are to be processed. First, the positive ions are 
treated, and then the electrons are treated by changing the sign of the 
potentials. There is a number density option and a current density option. 
The orbit quadrature loops are then set up, in which the energy and two 
velocity-angles are the variables of integration (summation). One option 
that of a single orbit; another option is that of a single energy (mon- 
energetic) in which only the angle sums are carried out. The quadratures 
are slightly different for the current and number density options. Each 
orbit is traced step-by-step backward in time until it terminates either 
on the metal surfaces or at “infinity" (outer boundary of the grid). At 
each step of the orbit, OENSTY calls INTERP and TRACK. It calls INTERP 
to find the potential and the components of the potential gradient (by 
interpolation within the grid) for use in the equations in motion. It 
calls TRACK to update the position and velocity components through the 
Newton equations of motion. In INTERP z and r are located by means of a 
table search, either in the North region or in the South region. In order 
to save time, after the first step in vihich the whole table is searched, 
the search is confined only to nearest-neighbor table entries. After the 
"box" containing z and r is identified, the potentials at the 4 corners 
of the box are used for a double-linear interpolation. 

TRACK, where the Newton equations of motion are used, monitors the 
"non-conservation" of energy, that is, the relative deviation between the 
assigned total energy E and the sum of the kinetic plus potential energies, 
where the latter two are subject to numerical errors as the orbit is propa- 
gated. The degree of energy non-conservation is kept under control by 
means of a time-step control governed by an input variable (STEP), Since 
the time-step control logic is probably not self-evident, this logic is 
explained as follows. 

During each time-step, the acceleration a is assumed to be constant. 
Therefore, if a is approximately constant within each grid zone (i.e., 
between two adjacent grid points in r and in z), the error (and, therefore. 
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energy loss or gain) would be that incurred upon crossing a boundary 

between two grid zones during a time step. Let a.| and a 2 be the constant 

accelerations in the first and second grid zones, respectively. Then the 
absolute value of the energy loss or gain, 1 aE|, would be given by 

|aE| - 1^2 " 1 ^2 

where $2 is the path length in the second zone. Now we know a^ and not ^ 2 * 

However, we estimate ^2 “3]! to be of the same order as |a^|. Moreover, 

we estimate $2 to be the larger of the two quantities 

|v^|At or |a^|(At)^ 

where At is the time-step interval. Hence, |aE| is estimated as the larger 
of the two quantities 

|a^v^|At or a^(At)^ 


Thus, if |aE| = fE is considered as known, the appropriate At may be esti- 
mated as the smaller of the two quantities 




In the program, the fraction f is identified with the input variable STEP. 
There is an additional control on At, to avoid taking too many steps, 
namely, a minimum value At = aS/(v^( where AS = .01 somewhat arbitrarily, 
representing of the order of 100 steps per unit length. For no electric 
field. At is taken to be the minimum value, with STEP replacing AS. 
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C.3 FORTRAN Input Variables 


The following variables are read by DWAKE as input: 

Card No, 1 

DATE = any identification (alphanumeric, 80 columns). 

Card No, 2 (integers only, 5 columns each) 

NCOLSN = number of r-valggs on North (wake) face of disk. 

MCOLSE = number of r-values East of disk (between disk and grid 
boundary), 

NCOLSS = number of r-values on South (front) face of disk. 
NROWSN = number of z-values. North region (between disk wake 
surface and grid boundary), 

MROWSS = number of z-values. South region (between disk front 
surface and grid boundary). 

Card No, 3 (8 fields of 10 columns each, per card) 

RHONl = r-values on North face of disk. 

Card No, 4 (8 fields of 10 columns each, per card) 

RHOE = r-values East disk. 

Card No, 5 (8 fields of 10 columns each, per card) 

RHOSl = r-values on South face of disk. 

Card No, 6 (8 fields of 10 columns each, per card) 

ZN = z-values. North region. 

Card No. 7 (8 fields of 10 columns each, per card) 

ZS = z-values. South region. 

Card No, 8 (8 fields of 10 columns each, per card) 

PHI = potential values at the grid points on North and South 
faces of disk. 

Card No. 9 (7 fields of 10 columns each) 

DEBYE = Debye number. 
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ALPH = Iteration parameter a. 

RBOUNO « r*va1ye of East boundary of grid. (If this Is greater 
than unity, r-values East are proportionally scaled 
between unity and RBOUND.) 

ZNBOUNO « z-va1ue of North boundary of grid. (If this is greater 
than zero, z-values North are proportionally scaled 
between zero and ZNB0<JND.} 

ZSBOUNO « z-va1ue of South boundary of grid. (If this is less 
than zero, z-values South are proportionally scaled 
between zero and ZSBOUND.) 

RUAKE » radius in wake such that additional iterations are 
applied to grid points lying within RWAKE, if option 
(ITALL >1 on Card No. ^O) is chosen. 

ZFRONT = negative z-value of front face of thick disk. If this 

is not zero, the disk is treated as having thickness equal 
to -ZFRONT. (Program TDWAKE.) 


Card ht (integers only, 5 columns each) 


ITS number of iteration cycles allowed (beyond the initial 
one). If this is zero, only one iteration is performed. 

IT = initial iteration cycle number (zero for new problem; 

greater than zero if continuing from previous problem). 

NEWPHI = zero for charge-density option; greater than zero for 
ion-density option. 

MAME = nonzero only if additional accuracy (more trajectories) 
desired for points in wake near axis. 

ITALL = nonzero if it is desired to apply more iterations to 

points in wake than to other points; indicates that all 
points are to be computed once every "ITALL-th" itera- 
tion, If this is zero, all points are computed every 
iteration. 


Card No. 11 (if IT = 0) Normally applies to number densities. 

NPRINT = printout option regarding orbit details (NPRINT = 0 for 
minimum normal printout). (Column 1.) 
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M0=0»>0 Indicates one, o>* more than one, space point to be pro- 
cessed (normally greater than zero on Card No. 11). 

(Column 2.) 

MC=0,>0 Indicates whether charge- or current-density is to be 
calculated (normally equal to zero on Card No. 11). 

(Column 3.) 

MA = number of quadrature intervals M, for polar angle a, 

d 

MA=0 indicates single orbit. (Columns 4 and 5.) 

MB = number of quadrature intervals for azimuthal angle 8. 
(Columns 6-10.) 

ME = number of quadrature intervals M^ for energy E. ME=0 
indicates single energy. (Columns 11-15.) 

STEP = orbit step-size control, for control of accuracy. 

(Columns 16-20.) 

RSAVE = initial value of r for single or'jit or for single space 
point. (Columns 21-30.) 

ZSAVE = initial value of z for single orbit or for single space 
point. (Colun;ns 31-40.) 

ALPHA = initial value of polar angle a for single orbit. 

(Columns 41-50.) 

BETA = initial value of azimuthal angle B for single orbit. 
(Columns 51-60.) 

EE = value of energy E for single orbit or for monoenergetic 
velocity distribution. (Columns 61-70.) 

XMSAVE = drift Mach number for Maxwellian with drift. (Columns 
71-80.) 

Card No. 12 (if IT=0) Normally applies to current densities. 

Same as above, except that MD=0 and HC>0. 

Additional Data Cards for Continuation 

If IT (iteration cycle number) is greater than zero on Card No. 10, 
that is, if the calculation is a continuation, then an array (normally of 
potentials) punched from the previous calculation is read in, after Card 
No. 10 but before Cards No. 11 and 12. 
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C.4 FORTRAN Output Variables (Normally Generated by Program) 

t 

Output generated by DWAKE (main program) 

Geometric and other input data, potential arrays, and iteration cycle 
nundier. 

Output generated by PRINT (called by FIELD) 

In the followinn, CN, CS, CE, CW, and C are the non-vanishing coeffi- 
cients in the Poisson difference-equation matrix: (initial cycle only) 

INDX = index of unknown potential (equation number) 

INDXN = index of North-neighbor potential 
CN = coefficient of North neighbor in equation number given 
by INDX 

INDXS = index of South-neighbor potential 
CS = coefficient of South neighbor in equation number given 
by INDX 

INDXE = index of East-neighbor potential 
CE = coefficient of East neighbor in equation number given 
by INDX 

IfJDXW = index of West-neighbor potential 

CW = coefficient of West neighbor in equation number given 
by INDX 

C = coefficient of central unknown in equation associated 

with INDX 

V = volume of space associated with central unknown in 

equation associated with INDX 

Output generated by SEIDEL (called by FIELD) 

X = solution of Poisson problem 

a one-dimensional array of potentials 

Output generated by DENSTY 

PHIN = two-dimensional array of potentials, North region 

PHIS = two-dimensional array of potentials. South region 
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NUMBER s Builier of trajectories (orbits) for each energy. (Off 
the axis, this is given by 4 on the axis, this is 
given by since only one valL-e of g is necessary by 
syaaetry. ) 

If the nunber-density option is chosen (MT=0), the following output 
is given cn « single line: 

K, KSAVE, ZSAVE, PHISAV = point niffi4>er, r-coordinate, z-coordinate. 

and potential, of point where density is 
calculated. 

DcitSA, KNST, CO * ion nuiriier-density, electron nuwDer-dansity, and 
charge (net) density. 

If corputations are done at a (Troup of space points {M0>0, usually 
for all the grid points), the output Includes the ion or charge density 
summary tabulation: 

kSV, ZSV, CDSW = r-ccKjrdinate of point, z-coordinate of point, 
and ion or charge density at that point. 

To monitor which trajectory takes the most steps, the following out- 
put ii given on a single line: 

MOSTPS, NSAVE, KES, uES. KBS, JBS, KAS, JAS = the largest number of 

steps taken by a trajectory, the index of the 
associated space point, and the velocity- 
coordinate indices of the trajectory (in the 
quadrature su«i) taking the largest number of 
steps. 

In addition, the ion densities are printed out as two-dimensional 
arrays, BflORTH and DSOUTH. 

If the current-density option is chosen (MC>0), the following output 
is given: 

For every value of energy E ca'iculatec .y the quadrat'jre formula, a 
line is printed, consisting of 
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NOESC « iMirtter of trajectories escaping 

NUNBER « total iM«i>er of trajectCMrles at the point of Interest 

FRACT « the fraction escaping (NOESC/NUNBER) 

E « the energy E 

(^S * the partial density in the sumatlon associated Mith 

energy E 

The following output is also given, as a single line: 

RSAVE, ZSAVE, PHISAV = r-coordinate, z-coordinate, and potential, 

of point kdiere current density is calculated. 
PARTCL, DEMST = "ion" or "electron," and value cf current 

o?nsity. 
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The listing for the TOUAKE program is given in the following pages. 



o o o o o o 


OVF’^LAYCOISKUS*0,8) . - - 

********** modification for finite disk thickness 

PROGRAM TOHAKEIINPUT, OUTPUT, TAPE6D=INPUT‘,TAPE61=0UTPUT,PUNCHJ . 

SATELLITE HAKE PROBLEM CUNSYHHETRIC OISK..-.DISKJUS) 

MIXING OF POTENTIALS FATHER THAN DENSITIES 

OPTIONAL ITERATION BY SWEEPING CROUPS OF POINTS IUPSTRE4N TO OOHNSTREA!^ 
MODIFICATION TO TREAT PILLBOX OR THICK DISK C>ZFRONT=THICKNESSl 

COMMON JJN,IIN,JJS,IIS,NT0T,RNTf5Gl ,ZNT(53) •RSTfSOl ,ZST(S0) , 

.1 RZI 5:0,21,PHIMI2C,2Cl,PHISf2C,20) ,CDC Sa:) .IFIRST.H 

COMMON/FLO/NCOLSN,MCOLSE,NCOLSS,NROMSN,X<5iO> *NROWSS,DcaYE,OEBYE2, 

1 RHON1(5CI,RHOE(50) ,RHOSU53) ,ZNt50), ZSC53) ,PHI (501 ,NGAP,NDISX 
C0MM0N/0EN/MPRINT,H0,MC,HA,H3,NE,STEP,RSAVc,ZSAV£,ALPHA,8ETA,EE, 

1 XMSAVE .RADIUS . 

COMMON/DENl/IT,HAME,NEHPHI,ISAVE,NGR,NGROUPf 533) ,DSAVEC 503) 

**^*K**** MODIFICATION, FOR FINITE DISK_tHICKNESS : 

i , ZFRONT 

DIMENSION DATE 120 1, Xl( 500 

C 

C FOR DENSITY AND CURRENT CALCULATIONS 

C INTEGER INPUTS ARE NPRINT,M0,MC,HA,H3,ME, NPRINT=3,1, 2, 3 F3R PRINTING 
C NO INTERMEDIATE DATA, ESCAPING INDICES, FIRST AND LAST STEPS, _ANO _ALL 
C STEPS. H0=0,l FOR ONE SPACE POINT IRSAYE, ZSAYE) OR SEVERAL SPACE 
C POINTS IREAO IN R,Z VALUES). HC=0,1 FOR CHARGE DENSITY OR CURRENT OEN- 
C SITY. MA=E MEANS ONE TRAJECTORY CREAO RSAVE, ZSAVE, ALPHA, BETA, EE) . ME-0 
C MEANS ONE ENERGY (READ EE). OTHERWISE, MA,M3,ME, DENOTE THE NUMBER 0 
C ALPHA-INTERVALS, BETA-INTERVALS, AND ENERGY-INTERVALS. 

C 

L=60 

M=61 

IPUNCH=0 

IFUNCH=1 ^ 

ISAVE=Q 

DEBYE2=0. 

i REA0(L,999a) DATE 

9998 FORMAT!20AA) _ 

IF(EOFfD) 99,15 

15 HRITEC«,9999) DATE 

9999 FORMAT C31H1W4SYHHE\ RIC DISK FIELD PROBLEM ,20A4) 

READfL,ll l) N COLS) -NCOLSE, NCOLSS,NROWSN,NROHSS 

NRCHSE=1 

JJN=NC0LSN»NC0LSE _ _ _ _ _ 

IIN=NR0HSN+1 

J «S=NCOLSS>NCOLSE _ _ 

^rS=NROWSS+l 

NOTSK=NCQLSNfNCOLSS _ 

NGAP=N0ISK*1 

NT9T=JJN»NR0HSN»NC0LSE<-JJS*K'0HSS 

RE-0(L,222) (RHONKJ) ,J = 1,NC0LSN) REPRODUCIBELIT? OF THE 

READCL,222) <RH0E U) , J=1 ,NCOLSE) ORIGINAL PAGE ^ POOR 

READ CL, 2 22) (RHOSl (J) ,J = 1,NC0LSS) 

RE!iO(L,222) (ZNCI) vI=l,NROWSN) 

REaoa,222) (ZSCI) ,I=1,NR0WSS) 

REA0(L,222) (PHlC J) , J=1 ,NOISK) 

PHI (NGAP) =.5* (PHI CNCOLSN) ♦PHI (NCOLSN^D). 

RA0IUS=RH0N1(NC0LSN) 
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9WAKE 



»••«••»«»» mooificaticm for finite disk thickness 

R:l 40CL.222I 0E3YE,4LPH,930UN0,ZNa0UN0,ZS00UN0.RHAK£,ZFR0NT 
WRITE (H, 229 }0EaVE«ALPH*R90UN0,2N30UN0,ZSa0UN0«R<IAK£fZFR0NT 
REAOlL.lllI ITSiIT.NEHPHItNANE.ITAU 
MRITEIN,lldl ITStlTtNEMPHItNAHEtlTALL 

ITMAXsTTS*IT 

DO 17 NslvNTOr 
XfNIsXlfNIsa. ' . 

OSAVEfNIsO. 

COIN1=0. ■' ' “ ■ 

17 CONTINUE 

NTOTAL=NTOT 

IFtlSAVE.GT.Q) REAOfL«52}NTOT,NTOT,fOSAV£(N} •Nsl,NTOT> 
IFIIT.EQ.OI GO TO 21 

IFtOESYE.GT.Q.) REA0ILt52}IT,NT0TALt C X INI •N^l ,NTOTAU 
IFIOEBYE.GT.Q,.ANO.NTOTAL.NE.NTOT> GO TO 19 
GO TO 21 

■ 18 WRIT£fN,668l 

668 FORMAT f////lX,*ilHTR0U8LE - INCORRECT NUMBER OF POINTS READ! 
GO TO 99 
C 

21 IFINEHPHI.EQ.a.ANO.OEBYE.GT.O.I 

1 HRITE(N,230I OEBYE«ALPH, ITHAX.IT, f N» XINI »N=1 ,NTOT| 
IFINEHPHI.GT.Q.AMO.OEBYE.GT.O.I 
1 HRITEfH»23Cl)0EBYE,ALPH.ITMAX,IT,tN« XI Nl •N=l»NTOTI 
C 

C RESCALE GRID LINE POSITIONS TO FIT WITHIN INPUT BOUNDS. 

C 

RATIO= IRBOUNO-RAOIUS) /fRHOE CMCOLSE) -RAOIUSI 
DO 23 J=l,NCOLSE ■ 

IFIRATIO.GT.Q.l RHOE I JI=RAOIUS ♦ RATIO*«RHOEI JI-RAOIUS) 

23 CONTINUE _ 

RATIO=ZNBOUNO/ZNC 11 

00 24 I=1,NR0WSN 

IFtRATIO.GT.O.) ZNm=RATIO«ZNfI> 

24 CONTINUE ^ 

RATI0=ZSB0UN0/ZS(NR0HSS1 

DO 25 I=1,NR0WSS 

IFIRATIO.GT.Q.l ZS f II =RATIO*ZS ill 

25 CONTINUE - - 

C 

WRITE lM,112rNCOLSN,NCOLSE,NCOLSS,NROHSN,NROWSS 
HRITEIH,223I CJ.RHONIIJI , J=1,NC0LSNI 
WRITEIH,224I IJ,RHO£f Jl ,J=1,NC0LSEI 
WRITECM,225I IJ.RHOSICJI ,J=1,NC0LSSI 
HRITEIM.2261 II.ZNIII ,I=1,NR0HSNI 
HRITEIM,2281 II.ZSII) ,I=1,NR0WSSI 

HRITEIH,229) 1 J,PHI f J1 , J = 1,NGAP| “ “ 

110 F0RMATI1X,3CHN0. OF ITS, IT, NEHPHI, NAME =, 415, 5X, 

1 16HALL POINTS EVERY , I3,24H-TH ITERATION AFTER IT=2I 

111 F0RMATI16I5I 

112 F0RMATC//1X,I3,25H COLUMNS CR-VALUESI NORTH/ 

1 1X,I3.25H columns IR-VALUESI EAST / 

2 " 1X,I3,25H COLUMNS IR-VALUESI SOUTH/ ~ 

3 IX,I3,22H ROWS (Z-VALUES) NORTH/ 

4 1X,I3,22H rows CZ-VALUESI SOUTHI 

220 FORMATI1X,6HOEBYE=,F10.5,5X,5HALPH=,F10.5, 
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1 1CX,7HR:.0UND » F7.2»5X*8HZN30UNO-> F7 .2«5X>8HZSaOUNO=. F7 .2, 
««««*«««*» M0DTFIC*4TICN FOR FINITE DISK THICHNESS 

2 5X. 6MRW«E=, F->-,2/lX,7M2FRONT=, F7.2I 

222 FORHATIf 10.5> 

223 FOPH«'^l/<^l5H_R-yALUES.HORTH/II3,lPEl5.%ll_ 

22 A F0RN/r(/«14H R-VALUES EAST/ II3,1PE15.^I I 

225 FO*» ’Tf//15M R-VALUES SOUTH/CI3* 1PE15. All 

226 F0RMATI//15H Z-VALUES NORTH/(t3«lP£l5. All 

228 F0RH4TJ//15H Z-WALUES S0UTH/II3.1PE15. All . 

229 FORMAT(*/2iW POTENTIALS ON DISK SURFACES/ ( 13 > 1PE15. Al 1 

230 _FORi*»T(/26HC . CONVENTIONAL ITERATION* , 

1 33H POTENTIALS HUH 0E8YE NUH8ER* F10.5*5X,6HALPH s* 

2 F«.0 r.,5X,6HITMAX = ,IA,5X,3HITs,IA/ClX,I3*lPE15.AII 

230^. FORHATf/26HCEXP-IN-POISSON ITERATION** 

1 30H POTENTIALS WITH OEBYE NUMBER* F10.5*5X*6HALPH. =* 

2 Fl2.«;,5X,6HITHAX=,IA,5X,3HITs,IA/(lX»I3.lPEl5.All 

5050 FORMAT flH0,IA,l&H ORDER POTENTIALI ' - 

5060 FORMATC1HO,IA,21H ORDER DENSITIES I 

5070 POR»*ATC1HO,IA,22H ORDER CURRENT OENSITYI 

52 FORMAH2I5,1P7E10,3/C 8E1G.3II 


DO 2 J=l«NCOLSN 
2_RNTfJl=RH0NHJ)_ 
DO 3 J=l.NCOLSE 
JPN=J»NCOLSN 
3 RNT(JPN)=RHOE{ji 
DO A I=l*NROHSN 
A ZNT(II=ZNCI1 

_ZNTIIINI=0. 

DO 5 J=l,NCOLSS 

5 RSTIJI=RHOSllJI_ 
DO 6 J=l,NCOLSE 
JPS=J*NCOLSS 

6 RSTtJPSI=RHOECJI 

ZSTC1I=0._ 

DO 7 I=l,MRbMSS 


IPS=I+1 


7 ZSTIIPS1 = ZSCII 

WRITE fM*231) (J*RNTIJ)*J=1*JJN) 

WRITE IN, 2321 IJ,RST( Jl *J=1 ,JJSI 

_HRITE (M.233I II*ZNTf I)*I=1,IINI 

WRITE tM,23Al f I, ZSTC I) ,1=1 *IISl 

231 FORMAT f //1X.27HAUGMENTEO R-VALUES, 

232 FORMAT (//1X,27HAUGMENTED R-VALUES* 

233 FORMAT f //IX, 27HAUGMENTEO Z-VALUES, 

234 FORMAT C//lX,27HAUGHENTEO Z-VALUES* 

C_ OUTPUT RHO ANO Z ARRAYS 

DO 71 I=l,NROWSN 

DO 71 J=1,JJN 

JRZ=(I-ll*JJN+J 
RZtJRZ,ll=RNT*J) 

71 RZ(JRZ,2)=ZNT(I1 

DO 72 J=l,NCOLSE 

JRZ=tNPOWSN*JJNI* j 
RZ(JRZ,l)=RNT(NCOLSN*Jl 

72 RZ(JRZ,2I=ZNT(IIN1 
NE=NROWSN»JJN»NCOLSE 


N0RTH/(I3,1PE15.4I) 

S0UTH/CI3,1PE15.4)I 

N0RTH/(I3,1PE15.4II 

S0UTH/(I3,1PE15.4I) 




ClBnj'IY OF THE 
■ IS Poor 


- 72 


TDWAKE 



00 73 r-ltNROMSS _ _ 

00 73 ' 

JR2=Ne*II-ll*J.JS^J 
RZCJRZ.ll-RSTUI 
73 R2CJRZ«2l=ZSm 
^ ___ 

NFPPs|NTOT/3aOI «>1 

OQ 85 IP=1,MFPP- “ 

WRITE IN. 801 

88 FORMAT C1H1,6X«1HI,^X«6H RfIl«%X.4HZIir I 
00 85 l3l,60 
K1=I*3C0* ClP-il 
K2=K1«60 
K3:=K2«60 
K4=K3*60 

K5=K4^60 ' ' 

IF(K5.LE.NT0T}HRITEIN,8IK1»RZIK1,1).RZ(K1,2) ,K2,RZ (K2,l) ,RZ(K2,2>, 

1 K3«RZ(K3«l)»RZfK3,2) «X4,RZ CK4,1I •RZCK4,2) ,K5,RZ(K5» 1) ,RZ C K3»2 ) 

IF fK5.LE.NT0n GO fO 85 

IFfK4.LE.NT0TIHRlTECH,8IKl,RZfKl,l) »RZfKl>2) ,K2tRZIK2,l} ,RZ(K2,21 , 

1 K3«RZ<K3«1)«RZ(K3«21 tK4«RZ (K4»l} •RZ(K4«2) 

IF (K4.LE.NT0n GO TO 65 

IF(K3. LE.NTOn WRITE (H,8)Kl,RZ(Kl,litRZCKl ,2> ,K 2,RZ(K2,1I ,RZ(K2,2) , 

■ '1 K3,RZ(K3.1)tRZ(K3,21 “ 

IF CK3.LE.NT0n GO TO 85 

IF(K2.LE.NT0T)HRITE(M,8IK1.RZ(K1,1) ,RZ(Klt2) .K2,RZIK2,ll »RZ(K2,2) 

IF CK2.LE.NT0n GO TO 85 

IFCK1.LE.NT0TIHRITEIM,8IK1,RZIK1,1I,RZCK1,2> 

85 CONTINUE 

8 FORMAT C5(I8«F10;^3;F8.3n 
RZ(NTOT+1*1)=0. 

RZ(NTOT+1,2)=0. ^ 

NGRPS=1 

IF(ITALL.LE.l) ITALL=1 ' 

IFIITALL.LE.il GO TO 1201 

C DEFINE GROUPS 1. 2t ETC.» IN THE HAKE, IN ORDER OF AXIAL DISTANCE FROM 0 
ZGROUP=ZN(NROHSNl 
NGR=1 

DO 12 N=1,NT0T ; 

NREV=NT0T-N»1 

NGROUP(NREV) = 0 ^ 

IF(RZ(NREV,2) .LE.O..OR.NREV.LE.JJN) GO TO 12 
IF(RZ(NREV,2) .NE.ZGROUP) GO TO 13 “ ’ 

IFfRZfNREV,!) .GT.RHAKEl GO TO 12 

NGROUP(NREV)=NGR 

GO TO 12 

13 ZGR0UP=RZ(NREV,2) 

TEMPORARY JUMP TO FORCE NGR=1 FOR ALL HAKE PTS. 

JUMP=1 

IF(JUMP.EQ.l) GO TO 12 
NGR=NGR«-1 
12 CONTINUE 

NGRPS=NGR 

C 

1201 CONTINUE 

REA0(L,666)NPRINTl,M0l,MCl,MAl,MBl,MEl,STEPl,RSftVEl,ZSAVEl, 
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1 ALPHA1,3ETA1,EE1,XHSAVE1 
666 FORHATt3Il.I2«2I5,F5.3t6F10 .51 

REA0a«6&6)NPRINT2«H02>HC2.HA2tHB2,HE2tST£P2tRSAtf£2*ZSA^E2, _ 
1 ALPHA2,3ETA2«EE2«XNSAV£2 « 

Tfirstso 


DO POTENTIALS 


10 00 30 NG::1«NGRPS 

NGR=NG 

IF(IfIRSf.EQ.Q.0R.IT.LE.2.0R.ITALL.LE.1.0R.H00fIT.ITALLI .EQ.3) 

t NGR=0 

CALL 0VERLAYt6L0ISKUS, 1*0,0) 

IF(NGR.LE.l) MRITECN.5C50) IT 

NG0=1 

11 IFCNGO.EQ.l) G0_T_0_%5. ■ 

DO 20 N=1,NT0TAL 

IF(IFIRST.EQ.O) X1IN)= XfN) 

X1CN)=ALPH» XfN) + ll.-ALPH)« XHM) 

20 X(M)= XI CM) 

45 00 50 I=1,NR0HSN 

00 50 J=lyJJN 

INDX=J»(I-1)*JJM 

50 PHIN(I,J)=XfINOX) 

DO 51 J=1,NC0LSN 

PHINCIIN,J)=PHIIJ) _ _ __ 

IFfIT.FQ.O) PHINCIIN*J)=0. 

51 CONTINUE 

JPLUSN=NC0LSN*1 

00 53 J=JPLUSN,JJN _ 

IN0X=IN0X»1 

53 PHINCIIN,J)=XCINOX) _ 

INOX=INOX-NCOLS£ 

JPLUSS=NC0LSS*1 

DO 54 J=1,NC0LSS 

NSU9VR=NDISK-fJ-l) __ 

PHrSfl,J)=PHI(NSUBVR) 

IFCIT.EQ.O) PHIS(1,J)=0. L 

54 CONTINUE 

DO 55 J=JPLUSS,JJS 

INOX=INOXU 
55 PHISC1,J)=XCINDX) 

DO 56 1=2, IIS 
DO 56 J=1,JJS 
INnX=INOX*l 

PHIS(I,J)=XCIN0X) 

■“56 CONTINUE 

IFCNG0.EQ.2) GO TO 500 
NG0=2 

HRITE(M,120) ITtNGR 
120 F0RMATC///,1X,23HP0TENTIAL ARRAY - NORTH, 5X,4HIT = , I 3,3X , 5 HNGR =, 

1 13) 

WRITECM,20G4) IRNTCJ),J=1,JJN) 

DO 91 1=1, IIN 

WRITE CM, 123) I.ZNTII), (FHIN(I,J) ,J=1,JJN) 

91 CONTINUE 






IS 


V:. 
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WRITE (N,122l 

122 FORHftT f///,lX,35H POTENTIAL ARRAY -SOUTH* //> 

WRITEIH.20SA} IRSTfJ),J=l*JJS} _ 

2004 F0RHAi(/lX,2HR3,16F6.4/C/3X*16F0.4ir 

00 93 1=1, IIS 

WRITE CM, 1231 I • ZST (I) , I PH IS 1 1 , J) • J=1 ; J JS) 

93 CONTINUE 

123 FORHATI/5H LINE,I4,5X,2HZ=,F8.V/(/7F16 .81) 

C 

GO TO 11 ~ ' " 

500 CONTINUE 

IFdPUNCH.CT.OI PUNCH 52, IT,NT0TAL, (X { Nr,N=r,NT3 TAL) 

IFIIT.GT.ITHAXT GO TO 1 

FIRST DO DENSITIES _ _ _ 

NPRINT=NPRINT1 

H0=H01 “ ~~~ 

HC=MC1 

HA=MA1 ■ “ " 

HB=H81 

HE=ME1 " ' 

STEP=STEP1 

STOP IN OENSTY IF'STEP LE ZERO. 

RSA\/E=RSAYE1 

ZSAVE=ZSAVE1 ~ 

ALPHA=ALPHA1 
8ETA=8ETA1 ~ 

EE=E£l 

XMSAVE=XHSAVE1 ^ 

CALL OVERLAYI6LOISKUS,2,0,Ol 
IFCIFIRST.EQ.0.OR.NGR.EQ.01 GO TO 31 

0 CONTINUE 

1 CONTINUE ' 

IFIRST =IFIRST^1 

WRITECH,6641 NPRINT , HO.MC.M A,MB, HE , STEP.RSA YE , ZSA ^EVALPhA; BETA ,EE, 
1 XMSAYE 

664 F0RMAT(1X,22HNPRINT,M0,HC,MA,HB,ME-,6I4/ 

1 1X,37HSTEP,RSAV£,ZSAVE,ALPHA,BETA,EE,XHACH=,7F10.5) 

WRITEIM,5C60I IT 

THEN DO CURRENTS 

NPRINT=NPRINT2 ‘ ‘ ' ' 

MD=M02 

MC=MC2 

MA=HA2 

M8=H02 

ME=ME2 

STEP=STEP2 

C STOP IN OENSTY IF STEP LE ZERO. 

RSAVE=RSAVE2 

ZSAVE=ZSAVE2 

ALPHA=ALPHA2 — - 

BETA=3ETA2 

EE=EE2 

XHSAYE=XHSAYE2 
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CALL 0VERLAYI6L0ISKUS,2.Q»6HR£CALLl 

MRITEfH«664) NPRlNTtNOtNC»HA»KB«NE.STEP,RSAVEtZSAV£i ALPHA. BETA >EE, 

i XMSAVE . ... . 

IF(HC2.GT.0I HRITCIN.5070I H 

IT=ITtl : ^ 

60 TO 10 

99 STOP 

END 
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SUBROUTINE ROOTIA«BtX| - - 

COHNON JJN.IINtJJS.IIStNTOT •RNTI501 .ZNT 1531 tRSTISOl t ZSTI50 ) • 

1 RZf 5a0»2»»PHINI26.Z0>.PHISCZCt20J .CO! 5Q3I f IFIRST«N 

C FIND ROOT OF X«B*EXPIXlsA. BY NEKTON NETHOO. 

KPRA=0 


KPRB=0 

EPSsl.C-6 . 


XOLD^XsO. 

KNAXslCQO 

* 

DO 100 K^l.KHAX 
XOLOsX 


KPRsK 

F=X ♦ B*EXPIX» - A 


FP=1. ♦ B’EXPIXI 
ox=o. 


IFfFP.GT.O.) OXs-F/FP 
X=XOLO ♦ DX 



0£LTA=OX 

IFIABStX) .GT.l.E-8) OELTA=OX/X 

IFIKPRA.GT.U) WRITE(H,1QQ0) K, A,B,X,OX , OELTA«F,FP 


1000 FORHATC1X,22HK,A,8,X,OX,OELTA.F,FP=,I5,1P7E14.4» 

IFfASSfDELTAt.LT.EPS) 60 TO 200 

100 CONTINUE 

WRITE CM, 9999) KHAX 

9999 FORMAT!// ///IX, 9HM0RE THAN, IS, _ 

1 40H ITERATIONS IN ROOT. HENCE PROGRAM STOP.) 

STOP 

C 

_ 200 _ CONTINUE J 

POELTAslOO. ♦DELTA 

IFCKPRB.GT.C) write (M,200f.) EP'?i,X,POELTA,KPR 
2000 F0RHAT!1X,35HC0KVERGENCE IN ROOT HtfHlN EPSILON=, 1PE9. 1 . IH . ,10 X, 
1 3HX =,E12.4,7H WITHIN, ElO. 2, llH 5ERCENT IN,I4,12H ITERATIONS.) 
RETURN 

END 
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OVERLAY fOISKUS,ttOI 
PROGRAH FIELD 


C UNSYMMETRIC DISK FIELD 

e. SEIOEL..METHOD 

C 

c'"’ ^ 

CONNOR JJN,IIN.JJS»IISfNTOT,RNTI50}tZNT(53) »RST fSO) • ZSK50 ) • 

1 RZ( 5GQ«2)«PHIN(20«2C) «PHISC2C,201 tRHANOt 5C3 ) tIFIRST.H 

COMMON/FLO/MCOLSN,NCOLSE,NCOLSS.NROHSN,X<5aO),NROHSS*OEaYEfO£BYE2, 

1 RHONl(50)«RHQEfSC) •RHOSK50) •ZN(5a>« ZS(50) .PHI (5 0) «NGAP,NOISK 

CONNON/A/CN«CS«CEtCU,C«Vt INOX, JSAT»RHO«Z*INDXN( 50G ) ,1NDXS ( . 500 ) 
l,INOXE( 5CO),INOXH( 5CC), CONST! 500,6) 

COMHON/OENl/IT,HAMc,NEMPHI,ISAVE,NGR,NGROUPC 500) ,OSAVEC 500) 

*♦♦♦*♦♦**• MODIFICATION FOR FINITE DISK THICKNESS 

l_,_ZFROMT 

C 

C ASSUME ASYMPTOTIC HONOFOLE 

ALPHAF(RR,ZZ)=-ZZ/!ZZ*»2 ♦ RR**2) 

8£TAF(RR,ZZ) j=-RR/(ZZ**2 ♦ RR**2) 


C 

C NTOT 


222 


223 


C 

C 

c 

c 

* 45 ' 


RESTORED AFTER HODIFICATION IN OENSTY 

NTOT=JJN*NROHSN ♦ NCOLSE ♦ JJS*NROHSS 

IF(IFIRST.EQ.O) GO TO 45 _ __ 

IFCNEHPHI.EQ.O.ANO.OEBYE.GT.a.) HRITE(M,222) OE8Y£,IT,NGR, 
CN,RHANOCN),N=l,NTOT) 

IFCNEWPHI. GT.O. AND. DEBYE. GT.O.) HRITE(M,223) DEBYE ,I T,NGR, 

(N,RHANO(N),N=l,NTOT) __ 

FORMATdHl/lSHOFIELD CALCULATION, lOX, 

41HINPUT=CHARGE DENSITIES WITH DEBYE NUM8ER = , FIO . 5, 

10X,4HIT =,I3,3X,5HNGR -.13/ 

(28X,I3,1PE15.4) ) 

FORMAKIHI/ISHOFIELO CALCULATION, lOX, 

41HINPUT= ION DENSITIES WITH DEBYE NUM9ER=.F1 0 .5, 

10X,4HIT =,I3,3X,5HNGR =,13/ 

C28X,I3,1PE15.4) ) _ _ 

NORTH ♦ NORTHEAST REGION 


FIRST POINT, FIRST LINE 


CONTINUE 


HEPil' )DUCIBILITY OF THE 
uRujINAL PAGE IS POQR ~ ~ 


JSAT=0 
1 = 1 

IFCIFIRST.EQ.O) HRITE(M,333) 

333 FORMAT!// ///25H NORTH ♦ NORTHEAST REGION///) 

IFdFIRST.EQ.O ) WRITE(M,334) I 

334 F0RMAT!//5H LINE,I3,93H ' N H 

1C E S) 

J=1 


INOX=J 

INDXN(INOX)=0 
lNOXS!INDX)=J>JJN 
INDXEiINOX)= INOX Vi 
INDXW(INOX)=0 
Z=ZN(T) 

RHO=RNT!J) 
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HN=0. 

HS= Z - ZNtI+1) 

HE=RNT J+D-RHO 
HW=0. 

ALPHA^ALPHAFCRHOt Z) 

CN=0. ■ 

CS= 0.125*HE**2/HS 

CE= HS/4. 

CVI=0. 

C- D.125*HE*«WC/HS ♦ 2.*HS/HE • ALPHA*HE) 

V= HS*HE**2/16. 

CAUL PRINT "5* 

> 

MIOOLE POINTS, FIRST LINE 

JHAX=JJN-1 ~ ~ 

00 5 J=2,JHAX 
■IN0X=J ■ 

INOXN(INOX)=0 

IN9XS(IN0X)=J^JJN 
INOXEtINDX) =IN0X«>1 
INDXH(IN0XI=IN0X-1 
Z=ZN(I) 

RHO=RNTfJ) 

HN=0. 

HS=Z“ZN(IH) 

HE=RNT(J + 1)“RH0 
HW=RH0-RNT(J-1) 

ALPHA=ALPHAF(RHO,Z) 

“ CN=0 . “ 

CS--0.5*CHE+HW»/HS*IRHO»CHE-HH)/4.) 

CE=0.5*HS/HE*CPHO+H£/2.) 

CW=0.5*H3/HH»(RHO-HW/2.) 

C='?.5*HS* (HE*HW)*(RH0/HE/HH+(1. -ALPHA* HS)/HS**2*(RH0> lH£-HW)/4 , 
V=0.25*HS*(HE*HW) *IRH0+CHE-HH)/4.) 

5 CALL PRINT 

LAST POINT FIRST LINE ' " 


J=JJN 

INOX=J 

INDXN(INDX)=0 

INDXS(INDX»=J*JJN 

INDXE(INDX)=0 

INDXW(IN0X)=IN0X-1 

Z=ZNtI) 

RHO=RNT(J) 

HN=0. 

HS=Z-ZNCI+1) 

HE=0. 

ALPHA=ALPHAF(RHO,Z) 

BETA=BETAF(RHO,Z) 

HW=RHO-RNT(J-l) 

CN = 0. “ 

CS=0.5*HW/HS* IRHO-HW/4.) 
CE = 0. 

CH=0.5*HS/HH*IRHO‘HW/2.) 
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C=3.5»nHM/MS+HS/MW-AtPHA*Hm* IRHO-HH/4 .) -HS* I8ETA*RHQ+ 0 .c’5) i 
V=0 . 25*HS^HW* f RHO-HH/4, I 

CALL PRINT . 

HIOOLE..LINES 

IHAX=NROMSN-l . .. 

IF(IHAX.LT.2) 60 TO 16 

DO 1C I=2,IMAX 

IF IIFIRST.EQ.Q) WRITE IN, 3341 I 

_ 00 10 J-ltJJN ' 

IN9XsJ*fI-l)*JJN 

INDXNCINOXl=INOX-JJM _ 

INOXS(INOX) sINOXf^JJN 

Z-ZNtll _ _ 

HN=ZNlI-ll-2 

HS=Z-ZN(m,l 

-RHO=RNTCJ) 

JG0=2 

IFfJ.EQ.l) JG0=1 ■ 

IFfJ.EQ.JJNJ JG0=3__ _ _ 

GO TO f6,7,d),JG0 

6 lNDXECINOX)=INOXtl 

INOXH(INOX)=0. 

HE=RNTCJ*1)-RH0 ^ _ _ _ 

HW=0. 

CALL LEFTCRHO,HN,HStHE,HH,CN,CS,CEfCH,C,V) 

GO TO 0 

_7 INOXEfINOX)=lNOX+l 

INOXH(INOXI=INOX-1 

HE-RNTfJ + lJ-RHO 

HH=RHO-RNTtJ-l» 

CALL MIDDLE (RHOtHN,HS,HE,HW«CN,CS,CE, CM, C»VJ 

GO TO 9 

8 1N9XE (INOXI=0 

INDXH(INOX)=INOX-i 

HE=0, __ 

HW=RH0-RNT(J-11 

3ETA=BETAFtRH0,ZI 

CALL RIGHT tRH0,HN,HS,HE,HM,3ETA,CN,CS,CE,CM,C,V) 

9 CALL PRINT 

10 CONTINUE 

LAST LINE ABOllE DISK SURFACE, NORTH-NORTHEAST REGION 

16 I = NROHSN ' “ 

IF (IFIRST.EQ.OI_HRITE (M,334l_I 

00 11 J=1,NC0LSN 
INOX=Jf II-1)*JJN 
INOXN(INDX>=INOX-JJN 
INDXS(INDX)=0 
INDXECINDXI=INDX+1 

IFCJ.GT.DINOXWIINDX) =IN0X-1 

IFCJ.EQ.1)INDXH(INDX)=0 

Z=ZN(I1 

RHO^RNTCJ) 

HN=ZN(I-1)-Z 
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HS=Z 
HF RNTf 

hm5Rho*rntij-ii 
IFt..CQ.ll Mliat. 

V. tJ.EQ.ll CALL LEFTfRH0,HH,HS*HE«HH«CN,CS»CE9CU«CtVI 

IF(J.Gr.ilCALLHIOOLECftHO,HN«HS*He,HII,CN»CStCE.CM«C;VI 

JSATsJ 

call PRIHT ~ ~ 

11 CQNTINIfC 

C * -- -- 

C LAST LINE AMAT FROH DISK SURFACE* NORTH>NORTHEAST REGION 

C ' ■ 

JNlN=NCOLSNn 

00 15 J=JNIN,JJN 

INCX=J*(I-1I*JJN 
INOXNf INCXI =INOX- JJN 
IMOXS CINOXI =INOX»NCOLSE 
IFIJ.LT.J^J IGO=l 
IFfJ.EQ.JJNl IGO-2 
INOXHIINOXI=INOX-1 

Z=ZNII} ' 

RHO=RNTfJI 
HN=ZNII-iJ-2 
HS=Z 

HR-RHO-RHTf J-11 
CO TO 112*13} * IGO 

12 INDXEIINOXI=IMOXH _ 

HE=RNTIJ*1I-P.H0 - 

call MIDDLE fRHO*HN,HS*HE*HH*CN*CS*CE*CM*C*Vl 
GO TO 1% " 

13 IN0XEfIN0XI=8 

HE=0 ■ 

8ETA=8ETAFIRH0*Z) 

CALL RIGHT*RH0*HN*HS*HE*HW,8CTA*CN,CS«CE*CM*C*VI 
1^ CALL PRIN1 

15 CONTINUE 

C 

C EAST REGION 
C 

lFfIFIRST.EQ.tJ} HRITECM*4G4I 

444 F0RMATC/////12H EAST REGION///! 
r=i 

IFfIFIRST.EQ.O} MRITEtN*334! I 

00 26 J=1*NC0LSE 

INDX=J»NRCfSN*JJN 

INDXNf INOX) =IN0X-NC0LSE 

INDXSIINOX)=INOX«^JJS 
R >J=RHOECJ) 

Z=0. 

;gc=2 

Ir'fJ.Ea.llJGO=l 

IFfJ.EQ.NCOLSE) JG0=3 

HN-7N(NR0HSN} 

HS--ZSfl} 

GO TO f20*21*22)*JGO 



OOO , OOO I I ooooo 


20 INnxCfINOX|sINOX»t . 

INOXUfIKOXIsO 

JSAT=NGAP 

H£=«IH0EIJ*11 - RHO 

HH:=RNO>RNTCNCOLSNI 

GO TO 23 

21 XNOXEfINOX)sINQX»l 

INOXHItNOXIsINOX-i 

HE^RHOeUni - RHO 

HM-RI*0 - RHOCU-ll 

GO TO 23 

22 iNOXEdNOXIM 

IN0XM(IN0X)=I«0X-1 _ 

HE=0. 

HH=RHO - RHOEIJ'll _ _ _ 

BETA=8ETAFCRHC,2I 

GO TO .2% 

23 CALL NIOOLEfRHO«HN,HS,HE«HH,CN,CS.CE.CM,C.V} 

CALL PRINT 

GO TO 26 

24 CALL RIGHTfRHO«HNtHS,HE,HM,BETA,CN,CS,CEfCRtC>VI 
CALL PRINT 

26_CONTINUE 


SOUTH *■ SOUTHEAST REGION 


IFCIFIRST.EQ.O) MRITECM,555I 

555 FORMAT I/////25H SOUTH ♦ SOUTHEAST REGION/ //I 

00 41 I=1,NR0USS _ _ 

IFCIFIRST.EQ.OI HRITEfH,334l I 
DC 41 J^l.JJS 

INOX=J4^CI-1)*JJS»NCOLSE«NROHSN*JJN 

_ RH0=RSTIJ1 

Z=ZSfIl' 

1G0=2 _ 

JG0=2 

IFd.EQ.ll IGO=l „ „ 

IF(I.EQ.NROHSS) IG0=3 

_ IFfJ.EQ.ll JGO=l 

IFIJ.EQ.JJSI JGO=3 
GO TC t28,29«30)« IGO 

FIRST LINE 

28 IN0XNCIN0XJ=IN0X-JJS 
IFU.LE.NCOLSS) INOXNlINOXI=a 
IFf J.LF.NCOLSS) JSAT=N0ISK*1-J 
INOXStINDX)=lNOX>JJS 

HN=-Z 

Hs=z - zsfini 

GO TO^l 

MIOOLE LINES 

29 IN0XN(IM0X1=IN0X-JJS 
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INOXS (INOXI sINOX« JJS 
HN=2SII-AI - 2 

Hs=2 - zscmi 

GO TO 31 

LASfLiME ^ 

3D lNOXN(INOX|sINOX-JJS 
INOXS t INOXI sQ 
HN=ZSII-1I - 2 
HS=0, 

ALPHA=-ALPHAF ( RHO * 2 r 
31 GO TO f 34.35, 361, JGO 

FIRST POIN*' 

34 IN0XEIIN0X)=IN0X*1 
INDXMfINOX)=0 
HE=RSTIJ*1I-RH0 

HH=a. 

IFfIG0.EQ.3l GO TO 33 

CALL LEFTCRHO,HN,HS,HE,HM,CK,CS,CE,CH,C,VI 
GO TO 40 



NIOOLE POINTS 

35 INOXE (INOXI -INOX^l 
IN0XH(IN0X1=IN0X>1 
HE=RSTCJ*1I-RH0 

HH=RH0-RST(J-1I 

IF(1G0.EQ.3I GO TO 33 

CALL NIOOLE (RHO.HN.HS, HE, HM,CN,CS,CE,CH,C,Vr 
GO TO 40 

LAST POINT 


36 iNOXEflNDXMQ 

INOXH (INOXI =IK0X-1 

HE=0. 

HH=RH0-RST(J-1I 

8ETA=BETAF(RH0,ZI 
~ IF (IG0.EQ.31 GO TO ^3 



IITY OF THE 
‘ ^ ^mR 


CALL RIGHTf RHO.HN.HS, HE, HU. BETA, CN,CS,CE,CH,C, VI 
GO TO 40 

33 GO TO (37, 30, 391, JGO 


LAST LINE, FIRST POINT 


37 CN=fl.l25*HE**2/HN 
CS=0, 

CE=HN/4. 

CH=0, 

C=0.125*HE*(HE/HN ♦ 2.*HN/HE - ALPHA*HEI 

V=HN*HE»*2/16. 

GO rO 40 

LAST LINE, MIDDLE POINTS 
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FIELD 
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38 CN=Q,5«IHE»HMI/HN« IRHO« 

CS=0. 

CE=0.5*HN/HE*IRHO«KE/2.t 

CK=0 .5»HN/HR* IRHO-HM/2. 

C= e . 5*HM* IHE^HM) • I RHO/HE/HH ♦ f 1 . -ALPHA »HN1 • C RHO » I HF-HH I / 4 . 1 / HH 

121 . .. 

VsC.25Q«HN*IHE«HH)«{RHO 4 (HE-HHI/%.1 

GO TO 48 ... 

LAST_LIME._LAST_POINI 

C 

39 CN=0.5*HII/HM*IRHO-HM/4.| _ _ 

CS=8. 

CE=0. 

CH=0.5*MN/HM*!RHO-HH/2.1 

C=6.5*IIHH/HM4HN/HM-ALPHA*HMIfL*RHO-HH/4,^|j^Hh»l3ELTA?RHQ_i:_0,251l 

V= 0 . 25*HN *HH* ( RHO-H H/4 . 1 

40 CALL PRINT _ 

4t CCNTIWE 

IF(IFIRST.GT.O) CALL SEIOEL 

EMO 
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FIELD 



SU8R0UTINE LEFTfRH0,HN,HS.HEVHHiCN,CS«CEf6MtCt VI 
CN=0.125*HE**2/HN -- 

CS=CN*HM/MS 

CE=Q.29*IHN4^HS) ■ ' ■ ' ~ 

CMsQ. 

C^C •25MHN*HSI *11 •♦0‘.5*HE**2/HN/HS) 
V-HE**2*CHN>HS)/16. 

RETURN " ' 

END 


SUBROUTINE MIOOLE f RHO,HN,HS.HE.HH,CN»CS*Ct*CH,C, VI 

CN=0.5*IHE+»«»/HM*IRHO+IHE-HH)/l».l 

CS=CN*HN/HS 

CE=0 . 5 * IHN*HS I /HE • f RHO*ME/k 

CH=3.5*<HN*HSI/HH*(RH0-HV/2,# 

C=Q.5*<HN*HS1*CH£+HHMCRH0/HE/HW+CRH0*(HE“HHI/4.I/HN/HSI 

V= 0 . 25* CHN+HSJ ♦IHE+HHI * (RHO + ( Ht-HHI _/4 . 1 

RETURN 

END 



SUBROUTINE RIGHT IRH0,HN,HS,HE»HH,3ETA,CN,CSf CE,CH,C,VI _ 

CN=0.5*HH/HN*CRHO-HH/4.I 

CS=CN*HN/HS 

CE=0. 

C H=(J . 5» C HN+ HS » /HH • ( RHO-HH/2 . ) 

C = 0,5*fHN+HS)*lHH/HN/HS*tRH0-HH/4.l*CRH0-HH/2Vl/HH' 0£TA»RHOI 

V=n .25* CHN*HS» *HH* CRH0~HH/4 .1 

RETURN 

ENO 
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SUBROUTINE PRINT - 

COHNON JJN,irN,JJStIIS,NTOT,RKTI5C)«ZNT(50)«RSTf5a) fZSTfSO I t 
1 RZI 500«2) •PHINCZCtZC) ,PHISf2C,ZQl ,RHANO( 5Q3I tlFIRST.H 
COMHON/FLO/ NCOLSN * NCOLSE, NCOLSS , NROHSN , X (5 J 0 I , NR0HS3 , OE8YE tOE3yE2 , 

1 RHONlCSaJ.RHOEISO) ,RHOSl(50).ZN(53l ZSI5 0J ,PHI (5 31 .NGAP.NOISlC- 

CONNON/A#CN,CS«CE,CH,C,V, INGX, JSAT,RHO,Z, INOXNI 50J)tlNOXS( 500) 

l,INOXE( 5eO),INOXH( 500).CONSTI 503,61 _ 

CP=0. 

CONST (INOX, 11= CN 

IFUINOXNCINOXl.EQ.a ) . AND. (CN.NE.O .) ) CP=CN 

CONST(INOX,2l= CS 

IF((INOXS(INDX).EQ.0 ) .AND. (CS .NE .3 . ) ) CP=CS 

CONST(TNOX,3J= CE . . 

CON-T(INOX,<*)= CH 

IF( (INOXH(INOX) .EO.O I . AND. (CW.Nc .0 . )) CP=CH 

C ********** TEHPORARY HELMHOLIZ EQUATION 

IF(OE8YE2.GT.O.) C=C* V/(OE BY£2**2» 

C ********** TEMPORARY HELMHOLTZ EQUATION 

C0NSTfIM0X,5)= C _ 

C0NST(IM0X,6I= V 

IF (CP.GT,O..ANO.OE8YE.EQ.e.) RHANQ (INOX) = CP*PHI IJSAT) 

IF (CP.GT.O..ANO.DEBYE.GT." .) RHAND(INOX) = RhANO ( IN0X)»V/0E3Y E**2 

1 ♦ CP*PHlfJSAT) 

IF (CP.EQ.O..ANO.DE3YE.GT.O .) RHANO(iNOXI = RHANO ( INOX) »V/0E9YE**2 
IF (IFIRST.GT.G) GO TO 3 

WRTTE(M,1) INOX, INOXN(INOX), CONST ( iNOX, 1 I , INOXH(INOX), CONST( 
1INDX,4), IN0X,C0NST(INDX,5) , INOXE(INOX), CONST(INOX,3), INOXSdN^ 
20X), C0NST(1N0X,2) , CONST (INOX, 6) 

_ 1 FORMAK / 6H POINT, I«t, 3H/C( , 14, 2H) = , IPEIO . 4 , 3H/C C , I U, 24 ) = , 

1E13.4, 3H/ri,I4,2H) = ,E10.4, 3H/C( , 14 , 2H) = , El 0 .4 , 3H/C ( , 14, 2 H) = , 

2E1C.4, 5H/VOL=,E10.^J _ 

IF (CP.NE.C.) WRITE (H,2) J-^AT, CP 

2 F0RMAT(31H COEFFICIENT OF POTENTIAL NO, (,I3, 4H) IS,F10.5)_ 

3 RZCINOX, 1J=RH0 

RZ(IN0X,2)^Z • 

RETURN 

ENO 
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100 


SUBROUTINE SEIDEL 

CONHON JJN,IIN«JJS,IIS»NTOT,RNT(50)«ZNTf5ai «RSTf5Q) »ZST(50), 

1 RZI 5CQ,21tPHIN(ZCt201,PHIS(2G,2G),RHANOI SCO ) • IFIRST»N 
COHHON/rLO/NCOLSN«NC' LSE*NCOLSS,NROHSN,X(5JOI*NROWSS,OEt)YE,OE3YE2, 
1 RHOH1CSO»,RHO£I5CO .TMOSHSOI •ZNtSO)* ZSC501 ,PHI (51 > ,NGAP,NOISK 
C0HH0N/A/CN,CStCE»Cl4«C*V« IHDX. JSAT,RHOtZtINOXN( SU:)«INOXS( 502) 
ItINOXEI SCOlfINQXHl 50C),CONST( 5C0>6) 

COMHON/OENl/IT. NAME, NEHPHI.ISAVE. NCR, NGROUPC 5CJ),OSAYEI 500) 
HOOIFICATION FOR FINITE DISK THICKNESS 

1 , ZFRONT 

RAOIUS-RHONI(NCOLSN) 

0MEGA=l,9 
EPS= 0.00001 

ITMAX=2000 ~ ~ 

ITCOUN =0 

IPROLO=0 * ■ 


IGO=l 

IF CNEHPHI .GT* 0 . AND .OEBYF'. GT .0 . ) WRITE C M. 103 ) 
F0RMATC/1X.44HH0DIFIED POISSON PROBLEM TO INCLUDE EXPCPHD) 
IFCIT.GT.O.ANQ.IFIRST.GT.Oi CO TO 2 
DO 1 K =l,NTOT 

1 XCK)=0. 

2 ITCOUN =ITCOUN H 
OELTAH=0. 

DO 3 K=l.MTOT 

X1=XIKI ~ '■ 

SN=C0NSTfK,l)/C0NSTIK,5) 

SS=C0NSTfK,2)/C0NSTCK,5) 

SE=C0NSTCK,3)/C0NSTCK,5) 

“ SH=C0NST(K»L)/C0NST(K,5) 

SR=RHAN0IK)/C0NST(K,5) 

INDXNK=INOXNCK) 

INDXSK=INOXSfK) 

INOXEK=INOXECK) 

INDXWK=INOXHCK) 

~ AA=SR 


-M4JOWJCIBILIT¥-OF-5IiK- 
OfilGINAL PAGE IS POOR 


IFf INDXNK.GT.O) A A = AA^SN*XI i:«0XNK) 
IF(INOXSK.GT.O) AA=AA>SS*XC INDXSK) 

IFCINOXEK.GT. 0) AA=AA*SE*X( INOXEK) 
IFfINOXHK.GT.O) A A= AA+SH*X( INOXHK) 
IFCNEWPHI.EQ.a.OR.DEBYE.EQ.O.) GO TO 30 
C MODIFICATION TO INCLUDE EXPCPHD POISSON PROBLEM 
BB=CC*NSTCK,6)/C0NST(K,5)/0EBYE**2 
CALL ROOTCAA.BB.XX) 

X(K)=XX 
GO TO 35 
30 XCKI=AA 

35 CONTINUE 


SET PHI=0 AT ZN BOUNDARY 

irCK.LE.JJN.ANO.CIFIRST.GT.O.OR.NEWPHI.GT.G)) XCK)=0 . 
********** MODIFICATION FOR FINITE DISK THICKNESS 

C SET Phi TO BODY POTENTIAL AT ADDITIONAL BODY POINTS IF ANY. 

IFCRZCK.l) .LE.RADIUS. ANO.RZCK.2) . GE. ZFRONT. AND. RZ (K. 2) .LT..”;,) 
1 XCK)=PHIC1) 

C 
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XCKIsOHEGA*XIK)*(l.-OMeGAI*Xl ... 

0EI.TA=A3SIX(KI-X1I 

IFtXl.NE.a.l 0EL1 A=ABS<(X(KI-X1I/X1I 

IFfOELTA.GT.OELTAHI OELTAH^OELTA 

3 CONTINUE 

IFdTCOUN.GT.ITMAXl HRITE(M,UI ITHAX 

IF(ITCOUN.GT.ITHAX) GO TO 9 

11 FO9MAT<////10H HORE THAN, lA .lOHiTERATIONS) 

IPR=ITCOUN/500 

IF(IPR.LE.IPROLO) GO TO 15 

IP90L0=IP.R 

GO TO 10 

15 IF (0ELTAM.G1.EPS) GO TO 2 

ITERATION FINISHED 


9 IGO=2 

10 NFPP=INTOT/300» ♦ 1 

DO 51 IP=1,NFPP 

WRTTE(H,22I ITCOUN,EPS,OELTAH,OHEGA 

DO 51 1=1,60 _ 

K1=I ♦ 300*CIP-1I 

K2=K1„> 60 

K3=K2 ♦ 60 

K4=K3 ♦ 60 _ 

K5=K4 ♦ 60 

IFIK5.LE.NTOT>HRITEIM,3333) Xl.XfKll ,K2,XfK2) ,K3,X (K3I ,K4,X tK4» , 

1K5,XCK5) 

IF(K5.LE.NT0T) GO TO 51 

IF CK4.LE. NT OT) WRITE (M, 33331 Kl,Xf Kll ,K2 , X(K21 , K3,X CK3) ,K4,X (K4) 
IFCK4,LE.NT0T1 GO TO 51 

IF(K3.LE.NT0T) WRITE IN, 33331 Kl.XCKl) •K2,X(X2) ,K3,X(K3) 

IF(K3.LE.NTOT) GO TO 51 

IFCK2.LE.NT0T)HRITEtM,3333) Kl.XCKlI ,K2,X(K2) 

IFCK2.LE,NTOTI GO TO 51 

IF{K1.L£. NTOTI WRITE C H, 33331 Kl.XIKl) 

51 CONTINUE . 

3333 FORMAH5(I8,Fl6.8n 
GO TO fl5,4),IG0 

22 FORHAT(15H1SOLUTION AFTER, 16, 2X.25HITERATI0NS WITH TOLERANCE, F 12 .8 

1,8X,18HMAXIMUN DIFFERENCE, F 12 . 8 , 8X, 6H0MEGA= , F8 . 5) 

4 RETURN ' ■ ' 

END 


♦ 
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OVERLAY (OISKUS,2«OI . . . _ 

PROGRAM OENSTY 

COMMON JJN,IIN,JJS,IIS,NO»RN(50I,ZN(5C1 tRS(5ai «ZS(53)»RSV( SOQ). 

1 ZSVt SQ0),PHIN(2G,2u),PHIS{2C,2a),C0SV( 530 ) • IFIRST ,H 

_ COMMON/DEM/NPRINT,HD,MC.MA,M3,ME,STEPtRSAV£,ZSAV£tALPHA,BEIA,EEL, 

1 XNSAVEtRAOIUS 

COMMON/OENi/IT,HAHEtNENPHl> ISAVEi NGR.NGROUPI SCdl •OSAVEC 5G0I 

MODIFICATION FOR FINITE DISK THICKNESS 

1 • ZFRONT I _ _ ' 

COMHON/P/ INT, I GN, JG N, IGS , JG S 

COMMON/T/X,Y,Z,R,XOOT,YOOT,ZOOT, PHI, PH IR, PHIZ, ERGO, DERG, C1»C2,E 

DIMENSION AfZl ,DNORTH (20 ,2u ) ,OSOUTH (2C , 231 

DIMENSION PARTCLC2>,PART1(21,PART2(2),FAT£C2),£ND1(2),EN02 (2> 

DATA PARV1/6H ION ,6H /,PART2/6H ELECT, 6HR0N / 

DATA EHD1/6H ABS0R,6HBED /,EN02/6H ESCAP,6HES / 

PI = 3.1415926536 

A 1 11 =-l. /SORT (3, 1 

A(2)=-A(ll 

MSTEP=10000 

HSTEP=20000 

IPRINT=1 

IF(MC.GT,O.OR.MGR«GT.l I IPRINT=0 

_NTOT=ND 

NPTS=ND 

IINM=IIN-1 _ _ _ 

IISM=IIS-1 

N1=IINH*JJN ^ _ ____ 

N2=NT0T-IISH*JJS 

00 5 I=l,IIN REPP^OTTCreiUTY OF THE 

DO 5 J=1,JJN original PAGE IS POOR 

5 DNORTHII,J)=0. 

DO 6 1=1, IIS 

00 6 J=1,JJS _ _ _ __ _ 

6 OSOUTH(I,JI=0. 

5 3 CONTINUE _ 

C DO ONE CHARGE DENSITY ORCURRENT DENSITY, OR 00 ALL 

IF(MD.EQ.G) NPTS=1 __ 

IFIMD.EQ.CI RSV(1I=RSAV£ 

IFtMD.EQ.O) ZSVU)=ZSAVE 
IF(MC.EQ.O) WRITE CM, 6661 IT,NGR 

IFCMC.GT.O) HRITE(M,667> IT 

WRITE CM, 664 INPRINT, MO, MC,MA, MB, ME,STEP,RSAVE,ZSAV£, 

1 ALPHA, BETA, EE, XMSAVE 

666 F0°MAT(lhl/lQHCD£NSITIES,5X,4HIT = , 13 , 3X ,5HNGR =,I3) 

667 F0RMATC1H1>9HQCURREKVS,5X,4HIT =,131 

664 FORMAT! X ,22HNPRINT,M0,MC ,MA,MB,ME=,6I4/ 

1 1X,37HSTEP,RSAVE ,ZSAVP ,ALPHA,PET.* ,EE,XHACH=,7F10 .5/1 

C 

IFCNPRINT.EQ.O) WRITEIM,66CI NPRINT 
IF(MPRINT.EQ.l) WRITECM,6611 NPRINT 
IFCNPRINT.EQ. 21 ’Rnt(M,662l NPRINT 
IF INPRINT.EQ. 31 WRITE<M,663I NPRINT 

66C F0RMATC1X,8HNPRINT =,I2,38H ;EANS NO trAjECTORY PRINTING * 

661 F0RMAT(1X,8HNPRINT =,I2,35H INDICES 0. ESCAPING TRAJEC DRIES ONLY! 

662 FOPMATCIX ,2HN"RINT =, 1 2, 3&H =FIKST *■ LASi ‘'’’EPS OF EACH TRAJECTORYl 

663 FO=>MAT C1X,8HNPRINT =,IE,36H EVERY STEP OF ALL TRAJEC'iORI. I ) 
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HOSTPS=0 

C 00 95 LOOP ENOS AT END OF PROGRAM ' “ 

IFCND.GT.e.ANQ.HC.EQ.OI NPTS=ND»1 

60 DO 95 Nsl,NPTS 

RSA\/E=RStffNI 

2SAVE=ZSVIN) 

IFfN.LE.Nll GO TO 7 

IFfN.GT.N2i GO TO 6 

7 DO 1 1=1, IIN 

IFfZSAVE.EQ.2Nmi IN=I 

1 CONTINUE 

DO 2 J=lV JJN 

IFfRSAVE.EQ.RNfJ) 1 JN=J 

2 CONTINUE 

IFfN.LE.Nll GO TO 9 

8 00 3 1=1, IIS ^ 

IFfZSAVE.EQ.ZSflll IS=I 

3~ CONTINUE 

00 4 J=1,JJS 
IFfRSAVE.EQ.RSfJJ 1 JS=J 
CONTINUE 

CONTINUE 

IFfN.LE.NZ.ANO.IN.GT.u.ANO. JN.GT.Ol ONORTHfIN, JN) =0SAV£(N1 
IFIN.GT.MI.ANO.IS.GT.O.ANO.JS.GT.G) OSOUTHCIS, JS) =0SAVE(N) 
IFfNGR.EQ.Ol GO TO 15 

IFflT.GT.O.ANO.N.LT.NPTS.ANO.NGROUPfNl .NE.NGR) COSVf N)=OSAVEfNl 
IFfIT.GT.O.ANO.N.LT.NPTS.AND.NGROUPfNi .NE.NGR) GO TO 95 
15 CONTINUE 

********** NOOIFICATION FOR FINITE OISK THICKNESS 

IFIRSAVE.LE.RAOlUS.ANO.ZSAVE.GE.ZFRONT.ANO.rZSAVE.LT.a.1 COSV(N) = 0 . 
IFfRSAVE.LE.RAOIUS.ANO.ZSAVE.GE.ZFRONT.ANO.ZSAVc. LT.o.l GO TO 95 
IF(HC.EQ.C.ANO.ISAVE.EQ.O) DSAVEfNl=0. 

MASAV£=MA 

MBSAVE=M8 ’ ' 

MESAVE=ME 

STEPSV=STEP 

INCREA=0 

IFfMC.GT.C.OR.MAHE.EQ.C) GO TO 20 

INCREASE ACCURACY NEAR AXIS 

IFfRSAVE.LE.RN(2) .ANO.ZSAVE.GT.O.l MA=ME=16 

IF(RSAVE.LE.RN(2> . AND. ZSA VE .GT .0 . 1 STEP=.05 

' IFfRSAVE.LE.RNt21 .ANO.ZSAVE.GT.0.1 INCREA=1 

20 CONTINUE 

C FIRST He no THE IONS ' - - . 

SCALE=1. 

PARTCLfll=PARTlfll 

PARTCL(2I=PARTH21 

C RETURN FROM END OF MAIN FOR ELECTRONS 

237 IF ISCALE.GT.O..ANO.N.EQ.1.ANO.IPRINT.EQ.1) WRITEt M,3000) 

1 (RN(J) , J=1,JJN) 

3000 F0RMAT(///,1X,24H POTENTIAL ARRAY - NORTH//1X, 2HR= , 16F8. 4/ 

1 (/3X,16F8.411 
C 

IFfSCALEiGT.O.) XMACH=XMSAVE 

IFfSCALE.LT.O. ) XMACH=0. 

PHTHAX = 0. 

DO 11 1 = 1, IIN 
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DO 10 . _ 

10 PHlNtI,JI=SCALE*PHIN(I,J) 

IF fPHIHAX.LT.PHINdtlll PHIHAX s PHINCI.il : 

IF(SCALE.EQ.>1..0R.N.6T.1.0R.1PRINT.EQ.G) GO TO 11 
HRITEIM.333I .I.ZNtl) . CPHINO .Jl . J=1.JJNI 

11 CONTINUE 

333 F0RHATC/5M LINE,II».5X.2HZ-,Fa.^/(/7F16 .8)1 

C 

IFCSCALE.GT.O.,ANO.N.EQ.1.ANO,IPRINT.EQ.1I WRITEIH,40Cai 

1 (RSI JI.Jsl.JJSI 

^^030 .FORMAT I ////24H P.OTENTIAL ARRAT_=_S0UTH//1X,EHR=:,16E8 .A/ 

1 C/3X.16F6.<»n 

00 13 isl.IIS 

00 12 J=1,JJS _ 

12 PHISCI,JI=SCAI.E*PHISCI,JI ^ 

IF(SCALE.EQ.-1,,OR,N.GT.1.0R.IPRINT,EQ.OLJLO.XQ_JJ 

HRITECM.333I I ,ZS I (I * (PHISC I.J) , J=1 . JJSI 

13 CONTINUE _____ 

C 

C SET UP SUMS OVER TRAJECTORIES 

C 

IF CMA.EQ.OI_GO T0_32 

JAMAX=2 

»BMAX=2 _ 

KAMAX=MA 

K8MAX=M8 ^ _ __ 

NUM8ER=MA*H8»4 

_C DC ONLY ONE BETA ON AXIS (SYMMETRY) 

IF (REAVE. EQ.O. ) K8MAX = l 

IF (RSAVE.EQ.n, ) NUMBER=MA*2 _ 

IF(RSAVE.EQ.O.) J8MAX=1 

IF (SCALE. G1 .0. .AND. N.EQ.l) WRITE(H,S6e> MA.MB.NUMBER 
668 F0RMAT(/1X,I6,16H ALPHA-IKTeRVALS,lX,I6. 15H BETA-INTERVALS, IX , I£ , 

_ 1_24H TRAJEC^’ORIES PER ENERGY) 

C 

IF(ME.EQ.C) GO TO 31 _ _ _ _ 

JEMAX=2 

KEMAX=HE _ 

IF (SCALE. GT.O. .AN. .N.EQ.l) WRITE(H,670) ME 

^570 F0RMATI1X,I6,47H E^ ERGY T^fERVALS, WITH 2 ENERGIES PER. INTERVAL//) 

C 

GO TO 33 

C 

C SINGLE ENERGY 
C 

31 -EMAX=1 

KEMAX=l'" 

IF(Sv.ALE.GT.C. .AN0.4.EQ.1) HRIT£(M,573) EE 
673 F0R’1AT(1X,31H MONOENERGETIC CASE WITH ENERGY, FlO .5//) 

C 

GO TO 33 

C 

C SINGIE TRAJECTORY ONLY 
C 

32 JAMAX=1 
JBMAX=1 
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JEH4X=l 

KAMAXn "* ^ 

KBHAXsl 

KEHAX=1 " ■■■; ■ *■ 

NUM8ER=1 

WRITE fM,669l ■AtPHAr'BETA.EE 

669 FORMATC 18H SINGLE TRAJECTORY/ 7H ALPHA=,F20 .18. 8H DEGREES/ 

1 6H BETAS, F20 .18. 6H DEGREES/ 8H ENERGYs.FlO .61 
C ® 

ALPHAsALPHA*PI/180." ^ “ 

BETA=BETA*PI/180. 

WRITE! M.665I ALPHA78ETA 

665 FORMAT! IX, 3H0R,/ IX, 6HALPHA=,F11 .8, 8H RADIANS/ IX, 5HBETA=,F10 
1.8, 8H RADIANS) 

C 

SINA=SIN!ALPHA)“ ' ~ “ “ 

COSAsCOS! ALPHA) 


SUM OVER ALPHA, BETA, AND ENERGY 


33 OENST=0. 

DO IQGl KE=1,KENAX 
00 1001 JE=1,JEHAX~ 

OENS=0. 

NOESCsO 

DO 1000 K8=l, KBMAX 
DO 1000 JB=1, JBMAX” 

00 lOOC KA=1, KAMAX 
DO 1000 JA=1, JAMAX“ 

INITIAL POSITION 

RsRSAVE — 

Z=ZSAVE 
X=R 
Y=0^ 

INT=0 ■ 

CALL IMTERP 
INT=1 " 

PHISAV=PHI 

ASSUME BOLTZMANN FACTOR FOR ELECTRONS '! OVERRIDE)"' 
IF!ABS(PHD.GT.500.I GO TO 96 
IF(MC.EQ.C. AND. SCALE. LT.O.) OENST=EXP! -PHI) 
tF(MC.EQ.O.ANO, SCALE. LT.O.) GO TO 96 

IF(MC,EQ.O.ANO.ISAVE.GT.fl) DENST=OS AVE !N) 

IFfMC.EQ.O.ANO.ISAVE.GT.O) GO TO 96 


IFCSTEP.LE.O.) HRITE!M,111I 
U F0RMAT!/////1X,L3HST0P DUE TO STEP LE. ZERO ***** ***** *****} 
IFtSTEP.LE.O.) STOP 

INITIAL VELOCITY 



SPEEO=0. 

IFtME.NE.O) GO TO <♦! 
E=EE 
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IF(E.LT.PHI) NRITE(H,67ii} KE, J£* KBi JB, KA» JA, BETAl * ALPHAl, E ,PHI SAW 
IFfE.LT.PHIl GO TO lOQl 

GO TO 40 

41 CE=CAfJEl ♦ FL0ATf2*KE-l-MEn/FL0AT|MEI 

Es C1.«CE)/I1.-CE) ♦ AMAXIIPHI ,C.I 

IFlXMACil.GT.l.) E = XHACH**2* fl.4CE)/ll.-CE» ♦ AHAXICPHI, J .) 
IFIE.tT.O. I MRITECM,674) KE*JE»KB« JJ.KA.JA.BETAl ,ALPHA1,E,PHISAV 
IFIE.tT.O. I STOP 

C , 

40 SPEEO*SQRT«E - PHI> 

. . IFIMA.EQ.CI GO TO 39 _ 

CA =(AfJA) f FLOATI2*KA • 1 • HA) ) /FLOAT IHA) 

IFIMC.EQ.OI COSA=CA 

POMER=XMACH 

IFtMC.EQ.O.ANO.POHER.GT.l.l COS A=-l .*2, • ( Cl .♦CAl /2 . I •♦POWER 

IFCMC.EQ.CI SINA=SQRT(1. - C0SA**2> 

IFCMC.GT.O) SINA=SQRT(.5*C1. ♦ CA)) 

IF(HC.GT.O) COSA=SQRT(l. - SINA**2) 

CBETA=(AIJB)^ FLOAT(2*KB - I - HB) ) /FLOATCHB) _ 

BETA=PI*C1.+ CBETA)/2. 

C . 

39 XOOT=SPEEO*SINA»COS(BETA) 

_ Y00T=SPEE0*SINA*SIN(8ETA) 

ZDOT=SPEED*COSA 

C1=STEP*E 

C2=SQRT(C1) 

KSTEP=Q 

ERGHAX=0. 

__ ERGO=0. _ 

OERGMX=0. 

OERG=0. 

ALPHAl = ACOSfCOSA)*180./Pl 
BETAl = BETA*180./PI 
ALPHA =ALPHA1 

BETA=BETA1 _ 

ZOLO=Z 

rFCNPRiNT,NE.2.AN0.NPRINT.NE.3) GO TO 34 
C PRINT INITIAL CONDITIONS OF TRAJECTORY 

WRITE CM, 674) KE , J E , KB , J8, KA , J A , BET Al , ALPHA L , E, PHISAV 
674 F0RMAT11X,3(I3,I2) ,F17.6,F14.8,2X,1P2E11.3,2X,46H =KE,JE, K3,JB, 
1 KA,JA, BETAl, ALPHAl , E,PHI) 


WRITECM,659) 

659 F0RMATC13X,115HSTEPS X Y Z XOOT 

1 YOOT ZOOT ERGMAX OERGMX ) 

WRITE(M,8e8) KSTEP,X,Y,Z, XOOT, YDOT, ZOOT _ _ __ _ 

888 F0RMATC13X,I5,1P6E11.3) 

TAKE A STEP 


34 CALL TRACK 

kstep=kstep*i 

IF(NPRINT.EQ.3)WRITE(M,868) KSTEP » X , Y , Z , XOOT , YOOT , ZOOT 
ir( ABSCERGMAX) .LT.ABS (EROO) ) ERGMAX -ERCO 
IF CABS (OERGMX) .LT . ABS (DERG) ) 0ERGMX = 0£RG 
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IF(KST£P.LT,HSTEPI GO TO 35 
WRITE (M,999l MSTEP 
WRITE CM, 97) KSTEP,N,KE, JE.KB, jq.kA , M 
999 FORMATCIOH MORE THAN, IlC , 5HSTEPSI 
STOP 

'35 R=SQRTCX*X»Y*Y) • 

IF (R,LT. RADIUS. AND. SIGNCl, ,Z| , ME. -IGNIl. ,20LD)) GO TO 36 
**•♦•♦**** MODIFICATION FOR FINITE OIS< THICKNESS 

IFCR.LT.RAOIUS.AND.Z.GT.ZFRONT.ANO.Z.UT.a.l GO TO 36 
IF(R.GT,RN(JJN).OR.Z.GT.ZNai .oR.Z.LT.ZSCIIS)) GO TO 37 
ZOLDsZ 

CALL INTERP • : 

GO TO 34 . 

PARTICLE IS A8SORBEO 


36 CONTINUE 

IF CNPRTNT .NE.2.ANO.NPRINT.NE.3) GO TO 1002 
FATEC1)=END1C1) 

FATE(2)=EN01(2) 

GO TO 374 


PARTICLE ESCAPES 


37 IF(NPRINT.EQ.l) GO TO 372 

IF (NPRINT.NE.2.ANO.NPRINT.NE.3) GO TO 373 
FATE(1)=EN02(1> 

FATE(2)=EN02 2) ' 

GO TO 373 

' HR ITE ( M, 674 ) KE , JE , KB , J8 , KA iJ A , BET A 1 , A LPH A 1 , E , PHIS A V 
.~/o NOESC=NOESC*l 

IFCHE.EQ.'') GO TO 374 
C 

CSANGL=2D0T/SQRTCXD0T**2 + Y00T»*2 -»■ 200T*»2) 

XPON=-2,»XMACH*SQRT ;E)*CSANGL - E - XMACH**2 
IFIMC.EQ.C) C0EFF1=SPEE0/FL0AT(NUM3ER) 

IF (MC.EQ.O.ANO.POWER.GT.l.) 

1 CO£FFl=COEFFl*POWER•(Cl.^CA)/^.)*•CPOHER-i.) 

IF(ABS(XP0N).GT,5CC.) G^ TO 17k 

IFCHC.GT.C) C0EFF1=SPEE0**2/FLCAT(NUHBER| 

OAOO=COEFF1*EXP(XPON) 

D£MS=OENS ♦ DADO 

C 

374 TF(NPRINT.NE.2.ANO.NPRINT.NE.3) GO TO 10C2 

WRITE CM, 889) F ATE , KSTEP . X,Y ,Z ♦ XOOT , YOOT , ZOOT , ERGHAX, OERGMX 
889 F0=>MAT(1X,2A6,I5, 1P8E11.3) 

1002 CONTINUE " 

IF(MOSTPS.GE.KSTEP) GO TO lOOC 

KES-KE 

JES=JE 

K3S=KB 

J3S=J8 

KAS=KA ■ 

JAS=JA 

NSaVE=N 

MOSTPS=KSTEP 
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IQQQ CONTimiC 

ENO OF 4MGLE SUM 

FR4CT=FtO« INOESC ) /FLOAT fmjNBERl 

C 

. IFCMPMMT.6T.O.OR.CMO.EQ.O.AMD.MC.GT.OII . . _ 

1 WRITE CM, 671 1 NOESC tNUHBER* FRACT t £»OENS 

671 F0RNATC18M RATIO ESCAPING 15, dH OUT OF *15, 

1 15H OR A FRACTION ,F13.8, 13H AT ENERGY £s,F13.8,4X, 

2 6HC0ENS»,1PE13.4,1H1J 

C 

IFCHPRINT.EQ.ei GO TO 5555 _ 

IF(N€.NE.8.ANO.NC.EQ.CI WRITE! H, 675) 

IFf* .NE.6.ANO.HC.GT.O) WRITEfN,676) 

675 FOR'ITIIX, 66H0ENS IS THE SUH OF OAOO=SPEEO*EXP(XPOM) /NUMBER OVER g 

IALL rlRECTIONS//) I 

676 FORHAMIX, 67H0ENS IS THE SON OF 0ftD0=SPEE0*»2*EXPIXP0N) /NUMBER OV- 1 

lER A ICHlSPHrRE//) ■ 

C 

5555 IFCME.EQ.8) GO TO ISOI 

IFfHC.EQ.O) C0EFF2=4,/SQRTfPI)/fl. - CE) **2/FLOAT CHE) 

IFCW.GT.D) C0EFF2=2./C1. - CE) **2/FLOATCHE) 

IFfXNACH.GT.l.) C0EFF2=C0EFF2*XMACH»»2 

OENST=OENST ♦ C0EFF2*0ENS _ . 

C 

1001 continue _ _ _ _ _ 

C 

TFCKE.EQ.C.ANO.HC.EQ.C) OEMST=SPEED»FRACT 

IFtHC.EQ.C.ANO.HC.GT.O) 0ENST=SPEE0»*2*FRACT 

96 IFCMC.GT.C) URITECH,677) RSAVE,2SAVE,PHISAV,PARTCL,0ENST 

677 F0RHATI/6H AT R=,F13.8, /H AND Z=,F13.8, 13H, THE POTENTIAL IS=, 
IF13.8/1X.20H ANO THE NORMALIZED ,2A6,2CH CURRENT DENSITY IS=,1PE13 
2.4//) 

_C _ 

IFISCALE.LT.l..'ANO.MC.EQ.b) GO TO 91 

IFCSCALE.LT.1..AN0.MC.GT.0) GO TO 90 

SCALE=-1. 

PARTCLC1)=PART2(1) 

PARTCLI2)=PART2C2) 

OENSA=OENST 

■ IF ( MC . EQ . 0 . ANOiiS A i/E.EQ . 0 ) OSA VE f N ) =OE NS A 

IF(N.LE.N2.ANO.IN.GT.C.AND.JN.GT.O) ONORTHCIN, JN) =OENSA 
IFCN.GT.Nl.ANO.lS.GT.C.ANO.JS.GT.C) DSOUTHC IS, JS) =OENSA 
IFCHO.GT.O.ANO.HC.EQ.O.ANO.N.EQ.NPTS) ONORTH (IIN,l):OENSA 
GO TO 237 

C RETURN TO BEGINNING ^^TRAJECTORIES FOR ELECTRONS 

■■'c ■" 

C CONTINUE IF IONS ANO ELECTRONS COMPLETED _ 

91 CO=OENSA-DENST 
CDSVfN)=CO 

IFfHC.EQ.O. ANO. NEHPHI.EQ.Q) DSAVE(N)=CO 
C SAVE ION DENSITY ONLY IF EXP(PHI) IS TO BE INCLUOEO_IN POISSON_ SOLUTION 
IF(MC.EQ.G.ANO.NEHPHI.GT.O) C0SVCN)=DENSA 
PHISAV=-PHISAV 

IFfMC.EQ.C) HRITE(M,672) N,RSAVE,ZSAVE ,PHISAV,0ENSA,DENST,C0 
672 FOPMATCIX, SHAT N=,I4,9H R,Z,PHI=, 1P2E1Q .2,E12.4, 
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IFIINCREA.GT.OI HRITE CH.6721} N,NA,M8>H£,ST£P 
6721 FORHATflX,12HFOR POINT H-t I3.4X, 

1 33HHA«HB«HE, ANO STEP ARE CHANGEO TQ»3I5,4H ANOfFia;$t 

2 23H FOR INCREASED ACCURACY/! 

OC DO 92 Isl.IIN 
DO 92 * 

92 PHIN(I,J}sSCALE«PHIN(I,JI 

DO 93 1=1, IK “ 

00 93 J=1,JJS 

93 PHISCI,J)=SCALE*fPHIS«I',Ji» 
restore HA,MB,H£,STEP 

NA=HASAV£ ' " 

NB=HBSAYE 
ME=MESAVE 
STEP=STEPSV 
95 CONTINUE 

IFfHO.GT.QI MRITEIH,666) IT,NGR _ 

IFCMO.GT.OI MRITE(H,94) f N,FSV(N) «ZSV(NI •COSVfNl , N=1,NPTS) 

94 F0RHATt/4X,lHN,8X,lHR,12X,lHZ,12X,4H0/C0/flX,I4,3F13.5)l 

C TRAJECTORY WITH HOST STEPS. PRINT INOICES(N,ANO K ANO J INOICES) 

WRITE CH, 97) M0STPS,NSAVE*KES,JES,KBS,J8S,KAS,JAS 
97 FORHATdX, 15*14,' 3113,12) ,34H =NOSTST£PS,N, KE,JE, K3,J3, KA,JA ) 
IFtIPRINT.EQ.O) GO TO 99 

HP.ITE1M,2O01) 

IFIISAYE.GT.O) MRITECM,2a03) 

2001 F0RMAT(1H1/1X,21H0ENSITY ARRAY - NORTH“ ) 

HRTTEIK,2004) CRNf J) , J=l, JJN) 

2004 ■ FORMAT I/1X,2HR=,16F8. 4/ C/3X,16F8T4)1 ^ 

2003 FORNAT(1X,40HOENSITIES READ IN RATHER THAN CALCULATED //) 

DO 100 1=1, IIN -- 

WRITEfH,333) I,ZNCI) , (ONORTHf I,J) ,J=1, JJN) 

100 CONTINUE 

HRITE(H,20Q2) 

2002 FORMAT (////IX , 21H0ENSITYARRAY -^SOUTH/7) 

WRITE(N,2Q04) (RS ( J) , J=l, JJS) 

DO 101 1=1, IIS ■ 

WRITE(H,333) I,ZS(I), (OSOUTH(I,J) ,J=1,JJS) 

101 CONTINUE 

99 CONTINUE 

END 


RI^RdDtJCffllLllY OP THE 
diiiGiNAL PAGE IS POOR 
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SUBROUTINE INTERP 

COMMON JJN,IIN,JjS,IIS,N 0«RN1SC )«ZNC5Q) •RSC531 ,ZS t501.RSV( 5QZi, 

1 ZSV( 5Qa),PHlNIZ0,20l,PHlSC20,2QI,COSlft 5CQI •IFIRST.H 
C0HH0N/P/INT«I6N»JGN,IGS,JGS 

C0HM0N/T/X,Y.Z«R»X00T,Y00T»Z00T«PHI«PHIR,PHIZ>ER30,0FRG»C1 fC2»E 

IGN=J&NsIGS^GSn ■ ■ 

NCHsQ 

IG02 s • ' ■ 

IF(INT.NE.O) IG02 - IGO 
IGO-2 

IF IZ.GE.O.) IGOsi 
" "lFfIG0.NE.IG02riNT's 0 
GO TO C1»2),IG0 

NORTH Z 

ASSUMING ZNCIINI LESS THAN OR EQUAL TO Z LESS THAN OR E QUAL TO ZERO 

1 IFfZ.EQ.ZNIin IG=2 
IFCZ.EQ.ZNIlll GO TO 103 
IFIINT.NE.Q) GO TO 100 
CO 1C 1=2, IIN 
IG=IIN-I»2 

IF f Z. LT. ZNCIG-in“GO TO ‘103 

10 CONTINUE _____ 

ACCEPT IF ZNCIG) LESS THAN OR EQUAL TO Z LESS THAN ZN(IG-l). 

100 IFfZ.GE.ZNIIG-in GO TO 102 

IFIZ.CE-ZNIIGirGO TO“10% 

101 IG=IG*1 

IF(Z.LT.ZNflG)) GO TO 101 
GO TO 103 

102 IG=IG-1 ~ " 

IFfZ.GE.ZNCIG-ll) GO TO 102 

103 NCH=1 ^ 

104 CONTINUE _ _ _ . _ 

NORTH R 

AS SUMING RNCli LESS THAN OR EQUAL TO R LESS THAN OR EQUA L TO ^NCJJ) 

IFfR.EQ.RNCJJN)) JG=JJN-1 
IF(R.EQ.RN(JJN)1G0 TO 153 
IF(INT.NE.O) GO TO 150 

DO 15 J=2,JJN - - - 

JG=J-1 

IF(R.LT.RNCJ))"GO~TO 153 

15 CONTINUE 

ACCEPT IF RNCJGJ LESS THAN OR EQUAL TO R LESS THAN RN(JGfl). 

150 IF (R.GE.RNfJGHn GO TO 152 

IF^R.GE.RNfJG)) GO TO 154 

151 JG=JG-1 

IFCR.LT,RN(JG»)GO TO 151 
GO TO 153 
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152 JG=JG»1 _ - _ 

IFIR.GE.RNtJG»ll)GO TO 152 

153 NCH=1 

154 CONTINUE 

C ^ 

C SET UP FRACTIONS 

C . 

OELZ=ZNfIG-l)-ZNIIGI 

OELR=RN«JG<-ll-RMCJGl 

FZ^IZ-ZNCIGII/OELZ 

FR=(R-RNCJGl)/OELR 

P22=PHINCIG-1,JG»1) 

P21=PHIN(IG-1,JG I 

P12=PHINIIG •JG^ll 

P11=PHINIIG ,JG > 

IGN=IG 

JGN=JG 

GO TO 500 

C 

C SOUTH Z 

C _ _ 

C ASSUMING ZSdlSl LESS THAN OR EQUAL TO Z LESS THAN ZERO. 

C _ 

2 IFflNT.NE.3l GO TO 200 

DO 20 1=2, IIS 

IG=IIS-I*2 

IFCZ.LT.ZSfIG-lJ ) GO TO 203 _ 

20 CONTINUE 


C ACCEPT IF ZSCIGI LESS THAN OR EQUAL TO Z LESS THAN ZSIIG-1). 

C 

IFfZ.GE.ZSfIGII GO TO 204 

201 IG=IG-H _ 

20C IFfZ.GE.ZSflG-inGO TO 202 

IFIZ.LT.ZSIIGMGO_TO_201 

GO TO 203 

202 IG=IG-1 

IF(Z.GE.ZS(IG-ll)GO TO 202 

203 NCH=1 _ __ _ 

204 CONTINUE ' ‘ ^ 

C _ ^ 

C SOUTH R 

C 

C ASSUMING RSCl) LESS THAN OR EQUAL TO R LESS THAN OR EQUAL TO RSCJJ 

C 

IFfR.EQ.RSfJJSn JG=JJS-1 

IF(R.EQ.RS(JJSn GO TO 253 

IF(INT.NE.O) GO TO 250 

DO 25 J=2,JJS 

JG=J-1 

IFCR.LT.RS(J)) GO TO 253 
25 CONTINUE 

C _ _ 

C ACCEPT IF RSiJGI LESS THAN OR EQUAL TO R LESS THAN RS(JG+1). 

C 

250 IFfR.GE.RSfJG+l) ) GO TO 252 
IFCR.GE.RSfJG) JGO TO 254 
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251 JG=JG-1 

IF(R.LT.RS(JG)}GO TO 251 
GO TO 253 
*252 JG=JG*1 

IF(R.GE.RSIJ6»l))GO TO 252 
' 253 NCH=1 “ ' • ■ 

254 CONTINUE 
C 

C SET UP FRaCTIONS 
C 

D€LZ=ZSIIG-ll-ZSaG) 

0£LR=RS(JG*1) -RSCJGI' ‘ 

FZ=(Z-ZSCIG))/OELZ 
FR=CR-9SfJGn/OELR 
P22=PHISCIG-1, JG*1) 

P21=PHISCIG-1,JG) 

P12=PHISCIG,4G*1) 

P11=PHIS(IG,JGI 

IGS=IG 

JGS=JG 

C 

C INTERPOLATE 
C 

500 IF(NCH.EQ.fll GO TO 501 
C 

C SKIP IF NO CHANGE IN PHI-BOX 

C 

D1=CP22-P12) /DELZ 
D2=(P21-P11) /OELZ 
D3=(P22-P21) /OELR 
D4=fP12-Pll) /OELR 

501 PHIZ=D2 ♦ FR* (01-02) 

PHIR=04 ♦ FZ» (03-04) 

PHI=P11 f FR*(P12-P11) ♦ FZ*(P21-P11) ♦ FR»FZ» (P22-P21-P12*P11) 
RETURN 

END ^ 
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SUaROUTIME TRACK 

CO»1HON/T/X,Y,Z.R.XOOT,YOOT,200T,PHI,PHIR,PHIZ*ER30,OERG,C1,C^,E 
0€RG=(PHI ♦ X00T*»2 + YOOT**2 ♦ ZDOT**2> /£”! • “ ERGO 
ERGQsERGO ♦ OERG 

VMAX=A9StX00T) ♦ ABS(YOOT) ♦ ABSIZDOT) 


STEP CONTROL 


IFIR.EQ.O.I PHIX=0. 

IFCR.EQ.fl,) PHIY=0. 

IF(R.EQ.O.) GO TO 1 
■ PHTX=PHIR*X/R ' 

PHIY=PHIR*Y/R 
1 SS = AMIM1(C2* Cl/</MAXl 

0T=SS/«A8SCPHIX1 ♦ A8S(PHIYI ♦ ABSfPHIZ) ♦ l.E-6» 
OT=AHAXl(OT, .01/VNAXI 

C THE FOLLOWING CARD IS FOR ZERO-POTENTIAL TESTS 
IF{PHIR,EQ.O..ANO.PHIZ.EQ.C.) OT=Cl/E/VMAX 
X=X+0T*fXC0T-,25*0T*PHIX» 
Y=Y«^0T»IYC0T-.25*0T*PHIY) 
Z=Z+DT*CZD0T-.25*0T»PHIZ) 

X00T=XD0T-,5*0T*PHIX 

YOOT=YOOT-.5*OT*PHIY 

ZDOT=ZCOT-.5*OT*PHIZ 

RETURN 

END 
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